Ask Your Question
0

plot_vector_field returning blank screen

asked 2015-08-03 14:41:37 -0500

jford1906 gravatar image

updated 2015-08-04 11:55:46 -0500

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 ...
(more)
edit retag flag offensive close merge delete

2 answers

Sort by ยป oldest newest most voted
0

answered 2015-08-05 16:21:24 -0500

jford1906 gravatar image

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)
edit flag offensive delete link more
0

answered 2015-08-04 06:04:42 -0500

This 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 ...
(more)
edit flag offensive delete link more

Comments

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.

jford1906 gravatar imagejford1906 ( 2015-08-04 09:19:19 -0500 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2015-08-03 14:41:37 -0500

Seen: 65 times

Last updated: Aug 05 '15