Ask Your Question

jford1906's profile - activity

2018-10-18 17:26:35 -0500 received badge  Popular Question (source)
2018-04-20 08:56:46 -0500 received badge  Popular Question (source)
2016-11-12 13:01:37 -0500 commented answer Parametric plot with piecewise input

Works, thanks!

2016-11-12 13:01:25 -0500 received badge  Supporter (source)
2016-11-04 03:01:38 -0500 received badge  Nice Question (source)
2016-11-03 14:33:56 -0500 asked a question Parametric plot with piecewise input

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

2015-11-22 16:34:58 -0500 asked a question Suggestions for improving speed

I'm working on visualizing a flow in R^3, and specifically focusing on how the flow behaves near a torus on inner radius 1 and outer radius 2, centered at the origin. The idea is to take a point (x,y,z), and look at this point using coordinates (r,z), in the plane containing (x,y,z) and the origin. Then I take equations in polar coordinates in this plane, and apply them to the point. Moving back into R^3, I take a rotation of the point, and add in the result of the flow restricted to the plane. Most of these conditions also require that I know the minimum distance from (x,y,z) to the torus mentioned.

This is taking a very long time, especially when I try to animate the flow. The last animation I tried took 20 hours, and that was before I found a flaw in the construction. This one works, but I'm hoping to optimize it before I try it again. Any suggestions on how to speed things up would be appreciated.

x,y,u,v = var('x y u v')

ep = 0.1 #step size
length = 350 #iterations

#initial conditions
omega_seeds = [(-1,-1,2),(1,-1,2),(-1,1,2),(1,1,2),(2,2,2),(-2,2,2),(2,-2,2),(-2,-2,2)]

twopi = 2 * RDF.pi()
c(r) = 2/3*(4**(r+1)-4)

def norm2((x,y)):
    """
    Norm of `(x,y)`
    """
    return float(sqrt((x^2 + y^2)))

def norm3((x,y,z)):
    """
    Norm squared of `(x,y,z)`
    """
    return x^2 + y^2 + z^2

def norm((x,y,z)):
    """
    Norm of `(x,y,z)`
    """
    return float(sqrt(norm3((x,y,z))))

def dist((x,y,z),(xx,yy,zz)):
    """
    Distance from `(x,y,z)` to `(xx,yy,zz)`
    """
    return norm((x - xx, y - yy,z-zz))

def atan2((x,y)): #generalized arctangent
    if x == 0 and y < 0:
        return 3*pi/2
    elif x == 0 and y > 0:
        return pi/2
    else:
        a = float(arctan(abs(y/x)))
        if x > 0:
            if y > 0:
                return a
            elif y < 0:
                return -a
            elif y == 0:
                return 0
        elif x < 0:
            if y > 0:
                return pi-a
            elif y < 0:
                return a-pi

#Functions to determine distances of a point from a specified torus
def Torus_dist((x,y,z),n):
    A = vector([0,c(n),0])
    if n%2 == 0:
        B = vector([x,y,0])
    else:
        B = vector([0,y,z])
    if A == B:
        return 4^n
    else:
        t = (2*4^n)/(norm(A - B))
        v = t*B+ (1-t)*A
        C = vector([v[0],v[1],v[2]])
    if dist(C,(x,y,z)) <= 4^n:
        return 0
    else:
        D = vector([x,y,z])
        t2 = 4^n/norm(C-D)
        v2 = t2*D + (1-t2)*C
        E = vector([v2[0],v2[1],v2[2]])
        return dist((x,y,z),E)

def dz(x,y): #vertical component of the flow restricted to the plane
    r = norm2((x,y))
    if r <=1 ...
(more)
2015-09-17 15:17:14 -0500 received badge  Self-Learner (source)
2015-09-17 15:17:14 -0500 received badge  Teacher (source)
2015-09-17 06:32:50 -0500 answered a question Restricting vector field to a subset of a rectangular region

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)
2015-09-14 15:01:05 -0500 asked a question 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)])
2015-08-08 06:57:41 -0500 asked a question Efficiently computing the distance from a point to a set

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.

2015-08-05 16:21:24 -0500 answered a question plot_vector_field returning blank screen

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)
2015-08-04 09:19:19 -0500 commented answer plot_vector_field returning blank screen

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.

2015-08-03 14:41:37 -0500 asked a question 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 ...
(more)
2015-06-15 12:34:59 -0500 asked a question Data type error when constructing a vector field

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'

2015-06-08 17:11:40 -0500 commented answer Unscaled arrows in a vector field

That's working, thanks!

2015-06-08 17:11:31 -0500 received badge  Scholar (source)
2015-06-08 16:31:35 -0500 received badge  Student (source)
2015-06-08 12:43:19 -0500 received badge  Editor (source)
2015-06-08 12:36:17 -0500 asked a question Unscaled arrows in a vector field

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