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.Thu, 17 Sep 2015 13:32:50 +0200Restricting vector field to a subset of a rectangular regionhttps://ask.sagemath.org/question/29440/restricting-vector-field-to-a-subset-of-a-rectangular-region/ I'm trying to restrict a vector field plot to a specific region of the plane. I found the code for the plot_vector_field3d function online, and just modified that accordingly, but I'm having trouble with the 2d case. plot_vector_field always picks points in a rectangular region, but there's a part of the region I want to skip. Here's what I have so far.
x,y = var('x y')
C = circle((0,0),1)
# 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()
def dy(x,y):
r = norm((x,y))
theta = arctan(y/x)
dr = cos(theta)*(1-1/r^2)
dtheta = -sin(theta)*(1+1/r^2)
if r<= 1:
return 0
else:
return dr*cos(theta)-r*sin(theta)*dtheta
def dx(x,y):
r = norm((x,y))
theta = arctan(y/x)
dr = cos(theta)*(1-1/r^2)
dtheta = -sin(theta)*(1+1/r^2)
if r<= 1:
return 0
else:
return dr*sin(theta) + r*cos(theta)*dtheta
VF = plot_vector_field((dx(x,y),dy(x,y)),(x,-2,2),(y,-2,2))
show(VF + C,aspect_ratio=1)
I'm not able to upload an image file, but if you run this, I'm trying to get rid of all of the vectors within the circle.
For comparison, here's the working code in 3d. This generates a flow on a torus, with the vectors only coming from points on or in the torus, instead of the usual plot_vector_field3d, which puts the points on an entire cube. dx, dy, and dz are the functions that give the correct vector components.
twopi = 2 * RDF.pi()
DT = [] # discretized tori
n = 20 # points in the boundary
m = 5 #layers in the torus
#Discretizes T0
for k in xrange(m):
for i in xrange(n):
DT.append((2*cos(i*twopi/(n)), 2*sin(i*twopi/(n)),0))
for j in xrange(n/2):
x0 = (2 + cos(j*twopi/(n/2)))*cos(i*twopi/(n))
y0 = (2 + cos(j*twopi/(n/2)))*sin(i*twopi/(n))
z0 = ((k+1)/m)*sin(j*twopi/(n/2))
DT.append((x0,y0,z0))
points = [vector(DT[i]) for i in range(len(DT))]
vectors = [vector((dx(*point), dy(*point), dz(*point))) for point in points]
max_len = max(v.norm() for v in vectors)
scaled_vectors = [v/max_len for v in vectors]
VF = sum([plot(v,color=hue(v.norm())).translate(p) for v,p in zip(scaled_vectors,points)])
Mon, 14 Sep 2015 22:01:05 +0200https://ask.sagemath.org/question/29440/restricting-vector-field-to-a-subset-of-a-rectangular-region/Answer by jford1906 for <p>I'm trying to restrict a vector field plot to a specific region of the plane. I found the code for the plot_vector_field3d function online, and just modified that accordingly, but I'm having trouble with the 2d case. plot_vector_field always picks points in a rectangular region, but there's a part of the region I want to skip. Here's what I have so far.</p>
<pre><code>x,y = var('x y')
C = circle((0,0),1)
# 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()
def dy(x,y):
r = norm((x,y))
theta = arctan(y/x)
dr = cos(theta)*(1-1/r^2)
dtheta = -sin(theta)*(1+1/r^2)
if r<= 1:
return 0
else:
return dr*cos(theta)-r*sin(theta)*dtheta
def dx(x,y):
r = norm((x,y))
theta = arctan(y/x)
dr = cos(theta)*(1-1/r^2)
dtheta = -sin(theta)*(1+1/r^2)
if r<= 1:
return 0
else:
return dr*sin(theta) + r*cos(theta)*dtheta
VF = plot_vector_field((dx(x,y),dy(x,y)),(x,-2,2),(y,-2,2))
show(VF + C,aspect_ratio=1)
</code></pre>
<p>I'm not able to upload an image file, but if you run this, I'm trying to get rid of all of the vectors within the circle.</p>
<p>For comparison, here's the working code in 3d. This generates a flow on a torus, with the vectors only coming from points on or in the torus, instead of the usual plot_vector_field3d, which puts the points on an entire cube. dx, dy, and dz are the functions that give the correct vector components.</p>
<pre><code>twopi = 2 * RDF.pi()
DT = [] # discretized tori
n = 20 # points in the boundary
m = 5 #layers in the torus
#Discretizes T0
for k in xrange(m):
for i in xrange(n):
DT.append((2*cos(i*twopi/(n)), 2*sin(i*twopi/(n)),0))
for j in xrange(n/2):
x0 = (2 + cos(j*twopi/(n/2)))*cos(i*twopi/(n))
y0 = (2 + cos(j*twopi/(n/2)))*sin(i*twopi/(n))
z0 = ((k+1)/m)*sin(j*twopi/(n/2))
DT.append((x0,y0,z0))
points = [vector(DT[i]) for i in range(len(DT))]
vectors = [vector((dx(*point), dy(*point), dz(*point))) for point in points]
max_len = max(v.norm() for v in vectors)
scaled_vectors = [v/max_len for v in vectors]
VF = sum([plot(v,color=hue(v.norm())).translate(p) for v,p in zip(scaled_vectors,points)])
</code></pre>
https://ask.sagemath.org/question/29440/restricting-vector-field-to-a-subset-of-a-rectangular-region/?answer=29484#post-id-29484Figured it out.
x,y = var('x y')
C = circle((0,0),1)
# 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()
def dy(x,y):
r = norm((x,y))
if r <= 1:
return 0
else:
if x == 0:
theta = 0
else:
theta = arctan(y/x)
dr = cos(theta)*(1-1/r^2)
dtheta = -sin(theta)*(1+1/r^2)
return dr*cos(theta)-r*sin(theta)*dtheta
def dx(x,y):
r = norm((x,y))
if r <= 1:
return 0
else:
if x == 0:
theta = 0
else:
theta = arctan(y/x)
dr = cos(theta)*(1-1/r^2)
dtheta = -sin(theta)*(1+1/r^2)
return dr*sin(theta) + r*cos(theta)*dtheta
ep = 0.3
def arrowplot(x,y):
n = norm((x - dx(x,y), y - dy(x,y)))*(1/ep)
if n == 0:
return ((x,y))
else:
return ((x+dx(x,y)/n), (y+dy(x,y)/n))
A = []
scale = 10
for i in range(scale+1):
x = -2 + i*4/scale
for j in range(scale+1):
y = -2 + j*4/scale
if (-2.1 <= arrowplot(x,y)[0] <= 2.1) and (-2.1 <= arrowplot(x,y)[1] <= 2.1):
a = arrow((x,y),arrowplot(x,y),arrowsize=1, color='black') #hue(norm((x - dx(x,y), y - dy(x,y)))*(1/ep)))
A.append(a)
show(sum(A)+C,aspect_ratio=1) Thu, 17 Sep 2015 13:32:50 +0200https://ask.sagemath.org/question/29440/restricting-vector-field-to-a-subset-of-a-rectangular-region/?answer=29484#post-id-29484