2021-10-24 20:50:08 +0100 received badge ● Famous Question (source) 2020-10-04 07:26:12 +0100 received badge ● Notable Question (source) 2020-07-20 20:46:25 +0100 received badge ● Notable Question (source) 2020-03-14 19:54:48 +0100 received badge ● Popular Question (source) 2020-01-07 12:50:17 +0100 received badge ● Popular Question (source) 2020-01-07 12:50:17 +0100 received badge ● Notable Question (source) 2018-10-19 00:26:35 +0100 received badge ● Popular Question (source) 2018-04-20 15:56:46 +0100 received badge ● Popular Question (source) 2016-11-12 20:01:37 +0100 commented answer Parametric plot with piecewise input Works, thanks! 2016-11-12 20:01:25 +0100 received badge ● Supporter (source) 2016-11-04 09:01:38 +0100 received badge ● Nice Question (source) 2016-11-03 20:33:56 +0100 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 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 ... 2015-09-17 22:17:14 +0100 received badge ● Self-Learner (source) 2015-09-17 22:17:14 +0100 received badge ● Teacher (source) 2015-09-17 13:32:50 +0100 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 22:01:05 +0100 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 13:57:41 +0100 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 23:21:24 +0100 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 16:19:19 +0100 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 21:41:37 +0100 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 ... 2015-06-15 19:34:59 +0100 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 "", line 1, in File "_sage_input_26.py", line 10, in 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 File "/private/var/folders/hs/vjm81xg9127bn9_6d1c4fp5w37f6j0/T/tmpTVrC6O/___code___.py", line 3, in 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 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-09 00:11:40 +0100 commented answer Unscaled arrows in a vector field That's working, thanks! 2015-06-09 00:11:31 +0100 received badge ● Scholar (source) 2015-06-08 23:31:35 +0100 received badge ● Student (source) 2015-06-08 19:43:19 +0100 received badge ● Editor (source) 2015-06-08 19:36:17 +0100 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