ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Wed, 05 Aug 2015 23:21:24 +0200plot_vector_field returning blank screenhttps://ask.sagemath.org/question/28751/plot_vector_field-returning-blank-screen/ I'm working on the level set at z=0 of a complicated vector field, but I only get a blank screen each time I try to run it. Any thoughts? The functions I'm calling for the vector field seem to return non-zero values on their own, but nothing appears when I ask for the plot.
This is actually a small part of a level set for a flow on all of R^3. In R^3, I take an infinite sequence of nested tori, and the differential equation at each point is described by an infinite series, with each term having a coefficient determined by the distance from that point to each torus in the space. c1 and c3 represent the inner and outer radii of the first torus. c4 represents the boundary of the second torus, which is 4 times as large as the first, and rotated from the xy-plane to the yz-plane, and completely contains the first torus. So that's why I'd discretized everything, since I'd already done that for the torus case, and was just trying to scale down. If I call the region between c1 and c3 T0, then h0(x,y) should return the minimum of 1 and 1-d((x,y),T0). If I call the union of the region inside of c1 and the region between c3 and c4 T1, then h1(x,y) should return the minimum of 1, d((x,y),T0) and 1-d((x,y),T1).
x,y,u = var('x y u')
def dist2(x,y):
return sqrt((x[0] - y[0])^2 + (x[1] - y[1])^2)
def setdist(x,S):
d = [dist2(x,S[i]) for i in range(len(S))]
return min(d)
C = parametric_plot((cos(u),sin(u)),(u,0,2*pi)) + parametric_plot((3*cos(u),3*sin(u)),(u,0,2*pi)) + parametric_plot((4*cos(u),4*sin(u)),(u,0,2*pi))
c1 = [(cos(i*2*pi/99*u), sin(i*2*pi/99*u)) for i in range(100)]
c3 = [(3*cos(i*2*pi/99*u), 3*sin(i*2*pi/99*u)) for i in range(100)]
c4 = [(4*cos(i*2*pi/99*u), 4*sin(i*2*pi/99*u)) for i in range(100)]
def h0(x,y):
if sqrt(x^2 + y^2)>=1 and sqrt(x^2 + y^2)<=3:
return 1
elif sqrt(x^2 + y^2) <1:
return 1-setdist((x,y),c1)
elif sqrt(x^2 + y^2) >=3 and sqrt(x^2 + y^2)<=4:
return 1-setdist((x,y),c3)
else:
return 0
def h1(x,y):
if sqrt(x^2 + y^2)>=1 and sqrt(x^2 + y^2)<=3:
return 0
elif sqrt(x^2 + y^2) <1:
return min(1,1-setdist((x,y),c1))
elif sqrt(x^2 + y^2) >=3 and sqrt(x^2 + y^2)<=4:
return min(1,1-setdist((x,y),c3))
else:
return 0
def dx(x,y):
return y*(h0(x,y))
def dy(x,y):
return (-1)*x*(h0(x,y)) + y*(h1(x,y))
VF = plot_vector_field((dx(x,y),dy(x,y)), (x, -4,4), (y, -4,4),plot_points=15)
VF + CMon, 03 Aug 2015 21:41:37 +0200https://ask.sagemath.org/question/28751/plot_vector_field-returning-blank-screen/Answer by slelievre for <div class="snippet"><p>I'm working on the level set at z=0 of a complicated vector field, but I only get a blank screen each time I try to run it. Any thoughts? The functions I'm calling for the vector field seem to return non-zero values on their own, but nothing appears when I ask for the plot.</p>
<p>This is actually a small part of a level set for a flow on all of R^3. In R^3, I take an infinite sequence of nested tori, and the differential equation at each point is described by an infinite series, with each term having a coefficient determined by the distance from that point to each torus in the space. c1 and c3 represent the inner and outer radii of the first torus. c4 represents the boundary of the second torus, which is 4 times as large as the first, and rotated from the xy-plane to the yz-plane, and completely contains the first torus. So that's why I'd discretized everything, since I'd already done that for the torus case, and was just trying to scale down. If I call the region between c1 and c3 T0, then h0(x,y) should return the minimum of 1 and 1-d((x,y),T0). If I call the union of the region inside of c1 and the region between c3 and c4 T1, then h1(x,y) should return the minimum of 1, d((x,y),T0) and 1-d((x,y),T1).</p>
<pre><code>x,y,u = var('x y u')
def dist2(x,y):
return sqrt((x[0] - y[0])^2 + (x[1] - y[1])^2)
def setdist(x,S):
d = [dist2(x,S[i]) for i in range(len(S))]
return min(d)
C = parametric_plot((cos(u),sin(u)),(u,0,2*pi)) + parametric_plot((3*cos(u),3*sin(u)),(u,0,2*pi)) + parametric_plot((4*cos(u),4*sin(u)),(u,0,2*pi))
c1 = [(cos(i*2*pi/99*u), sin(i*2*pi/99*u)) for i in range(100)]
c3 = [(3*cos(i*2*pi/99*u), 3*sin(i*2*pi/99*u)) for i in range(100)]
c4 = [(4*cos(i*2*pi/99*u), 4*sin(i*2*pi/99*u)) for i in range(100)]
def h0(x,y):
if sqrt(x^2 + y^2)>=1 and sqrt(x^2 + y^2)<=3:
return 1
elif sqrt(x^2 + y^2) <1:
return 1-setdist((x,y),c1)
elif sqrt(x^2 + y^2) >=3 and sqrt(x^2 + y^2)<=4:
return 1-setdist((x,y),c3)
else:
return 0
def h1(x,y):
if sqrt(x^2 + y^2)>=1 and sqrt(x^2 + y^2)<=3:
return 0
elif sqrt(x^2 + y^2) <1:
return min(1,1-setdist((x,y),c1))
elif sqrt(x^2 + y^2) >=3 and sqrt(x^2 + y^2)<=4:
return min(1,1-setdist((x,y),c3))
else:
return 0
def dx(x,y):
return y*(h0 ...</code></pre><span class="expander"> <a>(more)</a></span></div> https://ask.sagemath.org/question/28751/plot_vector_field-returning-blank-screen/?answer=28753#post-id-28753This is more a comment and code review than an answer, but since it's a bit long I'm using the answer field.
Below I'll use `r` for `sqrt(x^2+y^2)` and `C_R` for the circle of radius `R` centered at the origin.
1. In the definition of `C` you don't need parametric plots. Use Sage's built-in `circle` command.
2. It looks like `c1`, `c3`, `c4` are meant to be discretized circles of radius 1, 3, and 4.
In their definition, it seems that there should be no `*u` inside the `cos` and `sin`.
3. Are you using discretized circles on purpose, or as a workaround for continuous circles?
If it's a workaround for continuous circles, note that `dist((x,y),C_R)` is just `abs(R-r)`.
To see that, draw the circles `C_R` and the circle of radius `r` centered at `(x,y)`.
4. The quantity `d = setdist(x,S)` is always non-negative, so `min(1,1-d)` is always `1-d`.
So your functions `h0` and `h1` only differ if `1 <= r <= 3`, and in this case `h0 = 0` and `h1 = 1`.
Did you mean something else in the definition of `h1`?
Assuming discretized circles were a workaround for continuous circles, I would rewrite your code as follows.
# Use "real double field" floating-point values for efficiency
zero = RDF.zero()
one = RDF.one()
four = 4 * one
twopi = 2 * RDF.pi()
# Circles as graphic objects
C1 = circle((0,0), 1)
C3 = circle((0,0), 3)
C4 = circle((0,0), 4)
C = C1 + C3 + C4
# useful functions
def norm2((x,y)):
"""
Norm squared of `(x,y)`
"""
return x^2 + y^2
def norm((x,y)):
"""
Norm of `(x,y)`
"""
return norm2((x,y)).sqrt()
# h0 and h1
def h0(x,y):
"""
Return piecewise affine function of `r = norm(x,y)`
Let C_n be the circle of radius n, for n = 1, 3, 4.
With r = norm((x,y)) = sqrt(x^2+y^2), h0(x,y) is:
* r inside C_1
* 1 between C_1 and C_3
* 4 - r between C_3 and C_4
* 0 outside C_4
"""
r = norm((x,y))
if r < 1:
return r
if 1 <= r <= 3:
return one
if 3 <= r <= 4:
return four - r
return zero
def h1(x,y):
"""
Return piecewise affine function of `r = norm(x,y)`
Let C_n be the circle of radius n, for n = 1, 3, 4.
With r = norm((x,y)) = sqrt(x^2+y^2), h0(x,y) is:
* r inside C_1
* 0 between C_1 and C_3
* 4 - r between C_3 and C_4
* 0 outside C_4
"""
r = norm((x,y))
if r < 1:
return r
if 1 <= r <= 3:
return zero
if 3 <= r <= 4:
return four - r
return zero
def dx(x,y):
return y * h0(x,y)
def dy(x,y):
return -x * h0(x,y) + y * h1(x,y)
If discretized circles are relevant, varying the number of points may help illustrate your point, so use a parameter.
# discretized circles
n = 99
DC1 = [(cos(i*twopi/n), sin(i*twopi/n)) for i in xrange(n)]
DC3 = [(3*x, 3*y) for (x,y) in DC1]
DC4 = [(4*x, 4*y) for (x,y) in DC1]
We would also need to rewrite the functions `h0` and `h1` using these discretized circles.
For that, I would rewrite your auxiliary functions as follows:
# more useful functions
def dist2((x,y),(xx,yy)):
"""
Distance squared from `(x,y)` to `(xx,yy)`
"""
return norm2((x - xx, y - yy))
def dist((x,y),(xx,yy)):
"""
Distance squared from `(x,y)` to `(xx,yy)`
"""
return norm((x - xx, y - yy))
def setdist2(u,S):
"""
Distance squared from point to set
"""
return min(dist2(u,v) for v in S)
def setdist(u,S):
"""
Distance from point to set
"""
return setdist2(u,S).sqrt()
Once you answer some of the questions above, we can investigate the vector-field-plotting part.
x, y = var('x y')
VF = plot_vector_field((dx(x,y),dy(x,y)), (x,-4,4), (y,-4,4))
(C + VF).show(aspect_ratio=1)
Tue, 04 Aug 2015 13:04:42 +0200https://ask.sagemath.org/question/28751/plot_vector_field-returning-blank-screen/?answer=28753#post-id-28753Comment by jford1906 for <div class="snippet"><p>This is more a comment and code review than an answer, but since it's a bit long I'm using the answer field.</p>
<p>Below I'll use <code>r</code> for <code>sqrt(x^2+y^2)</code> and <code>C_R</code> for the circle of radius <code>R</code> centered at the origin.</p>
<ol>
<li><p>In the definition of <code>C</code> you don't need parametric plots. Use Sage's built-in <code>circle</code> command.</p></li>
<li><p>It looks like <code>c1</code>, <code>c3</code>, <code>c4</code> are meant to be discretized circles of radius 1, 3, and 4.</p>
<p>In their definition, it seems that there should be no <code>*u</code> inside the <code>cos</code> and <code>sin</code>.</p></li>
<li><p>Are you using discretized circles on purpose, or as a workaround for continuous circles?</p>
<p>If it's a workaround for continuous circles, note that <code>dist((x,y),C_R)</code> is just <code>abs(R-r)</code>.</p>
<p>To see that, draw the circles <code>C_R</code> and the circle of radius <code>r</code> centered at <code>(x,y)</code>.</p></li>
<li><p>The quantity <code>d = setdist(x,S)</code> is always non-negative, so <code>min(1,1-d)</code> is always <code>1-d</code>.</p>
<p>So your functions <code>h0</code> and <code>h1</code> only differ if <code>1 <= r <= 3</code>, and in this case <code>h0 = 0</code> and <code>h1 = 1</code>.</p>
<p>Did you mean something else in the definition of <code>h1</code>?</p></li>
</ol>
<p>Assuming discretized circles were a workaround for continuous circles, I would rewrite your code as follows.</p>
<pre><code># Use "real double field" floating-point values for efficiency
zero = RDF.zero()
one = RDF.one()
four = 4 * one
twopi = 2 * RDF.pi()
# Circles as graphic objects
C1 = circle((0,0), 1)
C3 = circle((0,0), 3)
C4 = circle((0,0), 4)
C = C1 + C3 + C4
# useful functions
def norm2((x,y)):
"""
Norm squared of `(x,y)`
"""
return x^2 + y^2
def norm((x,y)):
"""
Norm of `(x,y)`
"""
return norm2((x,y)).sqrt()
# h0 and h1
def h0(x,y):
"""
Return piecewise affine function of `r = norm(x,y)`
Let C_n be the circle of radius n, for n = 1, 3, 4.
With r = norm((x,y)) = sqrt(x^2+y^2), h0(x,y) is:
* r inside C_1
* 1 between C_1 and C_3
* 4 - r between C_3 and C_4
* 0 outside C_4
"""
r = norm((x,y))
if r < 1:
return r
if 1 <= r <= 3:
return one
if 3 <= r <= 4:
return four - r
return zero
def h1(x,y):
"""
Return piecewise affine function of `r = norm(x,y)`
Let C_n be the circle of radius n, for n = 1, 3, 4.
With r = norm((x,y)) = sqrt(x^2+y^2), h0(x,y) is:
* r inside C_1
* 0 between C_1 and C_3
* 4 - r between C_3 and C_4
* 0 outside C_4
"""
r = norm((x,y))
if r < 1:
return r
if 1 <= r <= 3:
return zero
if 3 <= r <= 4:
return four - r
return zero
def dx(x,y):
return y * h0(x,y)
def dy(x,y):
return -x * h0(x,y) + y * h1(x,y)
</code></pre>
<p>If discretized circles are relevant, varying the number of points may help illustrate your point, so use a parameter.</p>
<pre><code># discretized circles
n = 99
DC1 = [(cos(i*twopi/n), sin(i ...</code></pre><span class="expander"> <a>(more)</a></span></div>https://ask.sagemath.org/question/28751/plot_vector_field-returning-blank-screen/?comment=28757#post-id-28757I've added revisions above to clarify. I think I can get away without discretizing the circles for now, but I'll need to go back to it when I work out the 3d case. The point of getting a picture of this level set is to try and figure out the alpha and omega limit sets of the various trajectories in the 3d case.Tue, 04 Aug 2015 16:19:19 +0200https://ask.sagemath.org/question/28751/plot_vector_field-returning-blank-screen/?comment=28757#post-id-28757Answer by jford1906 for <div class="snippet"><p>I'm working on the level set at z=0 of a complicated vector field, but I only get a blank screen each time I try to run it. Any thoughts? The functions I'm calling for the vector field seem to return non-zero values on their own, but nothing appears when I ask for the plot.</p>
<p>This is actually a small part of a level set for a flow on all of R^3. In R^3, I take an infinite sequence of nested tori, and the differential equation at each point is described by an infinite series, with each term having a coefficient determined by the distance from that point to each torus in the space. c1 and c3 represent the inner and outer radii of the first torus. c4 represents the boundary of the second torus, which is 4 times as large as the first, and rotated from the xy-plane to the yz-plane, and completely contains the first torus. So that's why I'd discretized everything, since I'd already done that for the torus case, and was just trying to scale down. If I call the region between c1 and c3 T0, then h0(x,y) should return the minimum of 1 and 1-d((x,y),T0). If I call the union of the region inside of c1 and the region between c3 and c4 T1, then h1(x,y) should return the minimum of 1, d((x,y),T0) and 1-d((x,y),T1).</p>
<pre><code>x,y,u = var('x y u')
def dist2(x,y):
return sqrt((x[0] - y[0])^2 + (x[1] - y[1])^2)
def setdist(x,S):
d = [dist2(x,S[i]) for i in range(len(S))]
return min(d)
C = parametric_plot((cos(u),sin(u)),(u,0,2*pi)) + parametric_plot((3*cos(u),3*sin(u)),(u,0,2*pi)) + parametric_plot((4*cos(u),4*sin(u)),(u,0,2*pi))
c1 = [(cos(i*2*pi/99*u), sin(i*2*pi/99*u)) for i in range(100)]
c3 = [(3*cos(i*2*pi/99*u), 3*sin(i*2*pi/99*u)) for i in range(100)]
c4 = [(4*cos(i*2*pi/99*u), 4*sin(i*2*pi/99*u)) for i in range(100)]
def h0(x,y):
if sqrt(x^2 + y^2)>=1 and sqrt(x^2 + y^2)<=3:
return 1
elif sqrt(x^2 + y^2) <1:
return 1-setdist((x,y),c1)
elif sqrt(x^2 + y^2) >=3 and sqrt(x^2 + y^2)<=4:
return 1-setdist((x,y),c3)
else:
return 0
def h1(x,y):
if sqrt(x^2 + y^2)>=1 and sqrt(x^2 + y^2)<=3:
return 0
elif sqrt(x^2 + y^2) <1:
return min(1,1-setdist((x,y),c1))
elif sqrt(x^2 + y^2) >=3 and sqrt(x^2 + y^2)<=4:
return min(1,1-setdist((x,y),c3))
else:
return 0
def dx(x,y):
return y*(h0 ...</code></pre><span class="expander"> <a>(more)</a></span></div> https://ask.sagemath.org/question/28751/plot_vector_field-returning-blank-screen/?answer=28760#post-id-28760I figured something out, using a manual grid, and scaling the drawn arrows appropriately. I also modified it so any level set between z = -1 and z = 1 will work. Thank you for the help with cleaning up the code!
z = 0 #choose a z-value for the level set
x,y = var('x y')
# Use "real double field" floating-point values for efficiency
zero = RDF.zero()
one = RDF.one()
four = 4 * one
twopi = 2 * RDF.pi()
r1 = 2 + cos(arcsin(z)) #sets radii for the inner and outer circles
r2 = 2 + cos(arcsin(z) + pi)
# Circles as graphic objects
Cr1 = circle((0,0), r1,color='black')
Cr2 = circle((0,0), r2,color='black')
C4 = circle((0,0), 4,color='black')
C = Cr1 + Cr2 + C4
# useful functions
def norm2((x,y)):
"""
Norm squared of `(x,y)`
"""
return x^2 + y^2
def norm((x,y)):
"""
Norm of `(x,y)`
"""
return norm2((x,y)).sqrt()
# h0 and h1
def h0(x,y):
"""
Return piecewise affine function of `r = norm(x,y)`
With r = norm((x,y)) = sqrt(x^2+y^2), h0(x,y) is:
* r inside C_r1
* 1 between C_r1 and C_r2
* 4 - r between C_r2 and C_4
* 0 outside C_4
"""
r = norm((x,y))
if r < r1:
return r
elif r1 <= r <= r2:
return one
elif r2 <= r <= 4:
return four - r
else:
return zero
def h1(x,y):
"""
Return piecewise affine function of `r = norm(x,y)`
With r = norm((x,y)) = sqrt(x^2+y^2), h0(x,y) is:
* r inside C_r1
* 0 between C_r1 and C_r2
* 4 - r between C_r2 and C_4
* 0 outside C_4
"""
r = norm((x,y))
if r < r1:
return r
elif r1 <= r <= r2:
return zero
elif r2 <= r <= 4:
return four - r
else:
return zero
def dx(x,y):
return y * h0(x,y)
def dy(x,y):
return -x * h0(x,y) + -z * h1(x,y)
ep = 1.5
def arrowplot(x,y):
n = norm((x - dx(x,y), y - dy(x,y)))*(1/ep)
return ((x+dx(x,y)/n), (y+dy(x,y)/n))
A = []
scale = 15
for i in range(scale+1):
x = -4 + i*8/scale
for j in range(scale+1):
y = -4 + j*8/scale
a = arrow((x,y),arrowplot(x,y),arrowsize=3, color=hue(norm((x - dx(x,y), y - dy(x,y)))*(1/ep)))
A.append(a)
show(sum(A)+C,aspect_ratio=1)Wed, 05 Aug 2015 23:21:24 +0200https://ask.sagemath.org/question/28751/plot_vector_field-returning-blank-screen/?answer=28760#post-id-28760