| Hosted by CoCalc | Download
from sage.ext.fast_eval import fast_float var('x,y,z,a,b,c') f = y^4 + 5*x^4 - 6*x^2*y^2 + 6*x^3*z + 26*x^2*y*z + 10*x*y^2*z - 10*y^3*z - 32*x^2*z^2 - 40*x*y*z^2 + 24*y^2*z^2 + 32*x*z^3 - 16*y*z^3 ff = fast_float(f,'x','y','z') gradx = f.diff(x) grady = f.diff(y) gradz = f.diff(z) gx_ff = fast_float(gradx, 'x', 'y', 'z') gy_ff = fast_float(grady, 'x', 'y', 'z') gz_ff = fast_float(gradz, 'x', 'y', 'z') gmag = fast_float(sqrt(gradx^2+grady^2+gradz^2), 'x', 'y', 'z') pos_sphere = sqrt(16-x^2-y^2) pos_ff = fast_float(pos_sphere,'x','y')
(x, y, z, a, b, c)
%time N = 1000 R = 4 TOL = 1 # Carve up the xy plane into a (2N+1) x (2N+1) grid and try to find a z value of an intersectin at each lattice point our_points = [] for i in range(-N,N+1): if not i%100: print i x0 = R*float(i)/N x_used = False recent_y = None for j in range(-N,N+1): y0 = R*float(j)/N # Compute the z value on the sphere pos_z = pos_ff(x0,y0) # Compute the function value and the gradient. # If the gradient suggests we can get to an actual solution without changing by more than the size of our lattice, do it. pos_val = ff(x0,y0,pos_z) pos_grad = gmag(x0,y0,pos_z) pos_ratio = pos_val/pos_grad if abs(N*pos_ratio) < R/2: # Since we're not too far off, adjust our point to lie on our surface better and consider it done. our_points.append((x0-pos_ratio*gx_ff(x0,y0,pos_z)/pos_grad, y0-pos_ratio*gy_ff(x0,y0,pos_z)/pos_grad, pos_z-pos_ratio*gz_ff(x0,y0,pos_z)/pos_grad)) # Negatives now! neg_z = -pos_z neg_val = ff(x0,y0,neg_z) neg_grad = gmag(x0,y0,neg_z) neg_ratio = neg_val/neg_grad if abs(N*neg_ratio) < R/2: our_points.append((x0-neg_ratio*gx_ff(x0,y0,neg_z)/neg_grad, y0-neg_ratio*gy_ff(x0,y0,neg_z)/neg_grad, neg_z-neg_ratio*gz_ff(x0,y0,neg_z)/neg_grad)) print len(our_points)
-1000 -900 -800 -700 -600 -500 -400 -300 -200 -100 0 100 200 300 400 500 600 700 800 900 1000 7102 CPU time: 131.62 s, Wall time: 139.19 s
euclidean_distance = fast_float(sqrt((x-a)^2+(y-b)^2+(z-c)^2),'x','y','z','a','b','c') # There's probably something built in, but this'll do # Now we attempt to put the points in neighborly order MAX_DIST = 2 points_copy = our_points[:] neighbors = [[points_copy.pop(3000)]] # Start with not the first point in case the corner is weird counter = 0 while points_copy: if not counter%1000: print counter counter += 1 cur_point = neighbors[-1][-1] # Find the nearest point of the remaining ones nn = min(range(len(points_copy)), key=lambda i: euclidean_distance(*(points_copy[i]+cur_point))) dist = euclidean_distance(*(points_copy[nn]+cur_point)) if dist < MAX_DIST: # If it's not out of range, let's use it neighbors[-1].append(points_copy.pop(nn)) else: # Otherwise, check if it should loop... if euclidean_distance(*(neighbors[-1][-1] + neighbors[-1][0])) < MAX_DIST: neighbors[-1].append(neighbors[-1][0]) # ...and regardless start a new chain neighbors.append([points_copy.pop(nn)]) # Final loop check if euclidean_distance(*(neighbors[-1][-1] + neighbors[-1][0])) < MAX_DIST: neighbors[-1].append(neighbors[-1][0]) print len(neighbors), [len(block) for block in neighbors]
0 1000 2000 3000 4000 5000 6000 7000 2 [3552, 3552]
A = implicit_plot3d(f==0, (x,-4,4), (y,-4,4), (z,-4,4), plot_points=100, frame=false, opacity=.5) #Plot the surface C = implicit_plot3d(x^2+y^2+z^2==16, (x,-4,4), (y,-4,4), (z,-4,4), plot_points=80, frame=false, color="seashell", opacity=.2) the_visualization= sum((line(block, color="black") for block in neighbors), A+C) #Plot blocks of neighbors show(the_visualization.rotate((-1,-1,3),-5*pi/12).scale(2))
Error in lines 1-1 Traceback (most recent call last): File "/cocalc/lib/python3.9/site-packages/smc_sagews/sage_server.py", line 1230, in execute exec( File "", line 1, in <module> NameError: name 'f' is not defined