jford1906's profile - activity

Works, thanks!

I'm trying to plot a parametric function, which is defined piecewise. For some reason, the plot just jumps to the last piece of the parametric function and plots that.

Here's what I want to plot.

t = var('t')
r = 2

def f(x):
    if 0 <=x <=r:
        return (x, -r)
    elif r<x<=r + pi*r:
        return (cos(x-r), -(r + sin(x-r)))
    elif r + pi*r < x <= 3*r + pi*r:
        return (-x + 2*r + pi*r, r)
    elif 3*r + pi*r < x <= 3*r + 2*pi*r:
        return (cos(x - 3*r), -(r + sin(x-r)))
    else:
        return (x - 4*r - 2*pi*r, -r)

    parametric_plot(f(t), (t, 0, 4*r + 2*pi*r))

It just returns a straight line of length 4r + 2pi*r with y coordinate -r

Figured 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)

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)])

I'm working on a 3d vector field, where the derivative at each point is determined by a series, and each term of the series requires computing the distance from the point to a torus. Right now I'm only using 2-3 tori, and am discretizing each one, then computing a list of distances and taking the minimum. This is taking a really long time to run, and still isn't giving me the accuracy I need. Is there any way to improve efficiency in something like this?

The tori are defined via parametric equations, changing for each torus, and with moving centers.

Let c(n) = 2/3(4*(r+1)-4)

For every non-negative integer n, i,j in the range (0,2*pi), and k in {1/10, 2/10,...,9/10,1}, T_n is given by

x = 4^n*(2 + cos(j)*cos(i)
y = 4^n*(2 + cos(j))*sin(i) - c(n)
z = 4^n*k*sin(j)

if n is even and

x = 4^n*k*sin(j)
y = 4^n*(2 + cos(j))*sin(i) - c(n)
z = 4^n*(2 + cos(j)*cos(i)

if n is odd.

So far, I've been discretizing the tori using these equations, computing a list of Euclidean distances, and taking the minimum.

Here's the full equations for determining the vector field. First, I make a list of points for each torus. For any p = (x,y,z) in R^3, define a sequence of functions, with

h_0(p) = min(1, 1-d(p,T_0))
h_n(p) = min(1, d(p,T_(n-1)), 1-d(p,T_n))

Additionally, define

G_n(p) = (y,-x,0)

if n is even and

G_n(p) = (0,-z,y)

if n is odd.

Finally, for each point p in R^3, the derivative at p is given by

dp = sum(G_n(p - (0,c(n),0))*h_n(p))

over all non-negative integers n. So far I've been restricting myself to the n = 0,1,2 case, but even then, I generally can't get the field to display unless I make the discretization really large, and that doesn't show the detail I need to help me find out where the alpha and omega limit sets are.

I'm not having any problems writing the code that will make these calculations, I'm just wondering if there is a more efficient way, rather than taking a list of points for the torus, and just computing the minimum over and over again.

I 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)

I'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.

I'm trying to get a 3d vector field, using a smooth bump function, and I keep getting this error. I've searched around and can't figure out what the issue is.

x,y,z,r,theta = var('x,y,z,r,theta') 
ep = 0.1
f(x,y) = (x^2-1)*(y^2-1)*x*y
g(x,y) = (x^2+(y-3/8)^2-1/64)*(x^2+(y+3/8)^2-1/64)

def beta(t): #bump function on [0,ep]
    if (t>=0) and (t<=ep):
        return numerical_integral(e^(1/(x*(x-ep))),0,t))[0]/(numerical_integral(e^(1/(x*(x-ep))),0,ep))[0]
    else:
        return 0

def v(x,y): #returns 1 near singular points
    d = min((x^2 + (y-1/4)^2), (x^2 + (y-1/2)^2))
    if sqrt(d)>ep:
        return 0
    else:
        return abs(1-beta(sqrt(d)))

def Transform(x,y):
    r = sqrt(x^2 + y^2)
    return r

Y0_0 = f(x,y)/sqrt(f(x,y)^2+g(x,y)^2)
Y0_1 = g(x,y)/sqrt(f(x,y)^2+g(x,y)^2)

Y0 = plot_vector_field((Y0_0,Y0_1),(x,-1,1),(y,-1,1),plot_points=15) #returns a 0-annihiltor

ex(r,theta)=-r*sin(theta)
ey(r,theta)=r*cos(theta)

def Y1_0(theta, x, y):
    r = Transform(x,y)
    return -v(x,y)*ex(r,theta)

def Y1_1(theta, x, y):
    r = Transform(x,y)
    return -v(x,y)*ey(r,theta) + Y0_0(x,r)

def Y1_2(theta, x, y):
    r = Transform(x,y)
    return Y0_1(x,r)

Y1 = plot_vector_field3d((Y1_0, Y1_1, Y1_2),(x,0,2*pi),(y,-1,1),(z,-1,1),plot_points=15) #returns a 1-annihilator

Traceback (most recent call last): File "<stdin>", line 1, in <module> File "_sage_input_26.py", line 10, in <module> exec compile(u'open("___code___.py","w").write("# -- coding: utf-8 --\n" + _support_.preparse_worksheet_cell(base64.b64decode("cGxvdF92ZWN0b3JfZmllbGQzZCgoWTFfMCwgWTFfMSwgWTFfMiksKHgsMCwyKnBpKSwoeSwtMSwxKSwoeiwtMSwxKSxwbG90X3BvaW50cz0xNSk="),globals())+"\n"); execfile(os.path.abspath("___code___.py")) File "", line 1, in <module>

File "/private/var/folders/hs/vjm81xg9127bn9_6d1c4fp5w37f6j0/T/tmpTVrC6O/___code___.py", line 3, in <module> exec compile(u'plot_vector_field3d((Y1_0, Y1_1, Y1_2),(x,_sage_const_0 ,_sage_const_2 *pi),(y,-_sage_const_1 ,_sage_const_1 ),(z,-_sage_const_1 ,_sage_const_1 ),plot_points=_sage_const_15 ) File "", line 1, in <module>

File "/Applications/Sage-6.1.1.app/Contents/Resources/sage/local/lib/python2.7/site-packages/sage/plot/plot3d/plot_field3d.py", line 90, in plot_vector_field3d return sum([plot(v,color=cm(v.norm()),**kwds).translate(p) for v,p in zip(scaled_vectors, points)]) File "/Applications/Sage-6.1.1.app/Contents/Resources/sage/local/lib/python2.7/site-packages/matplotlib/colors.py", line 589, in __call__ lut.take(xa, axis=0, mode='clip', out=rgba)

TypeError: Cannot cast array data from dtype('O') to dtype('int64') according to the rule 'safe'

That's working, thanks!

I have a vector field I'm plotting, but the arrows are scaled so small near the singular points, it's unclear which are sinks, sources, and saddles. Is there a way to change the arrow scaling so they are all the same length and the behavior is clearer?

x,y = var('x,y')
VF = plot_vector_field(((x^2-1)*(y^2-1)*x*y, (x^2+(y-3/8)^2-1/64)*(x^2+(y+3/8)^2-1/64)), (x, -1,1), (y, -1,1))
points = point([(0,1/4),(0,1/2),(0,-1/4),(0,-1/2)],color='red', size=25)
VF + points