Instantly share code, notes, and snippets.

# dhermes/implicitizing_curves.pdf Last active Sep 17, 2018

Discussion of parametric curves; particularly classification and implicitization
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
 import matplotlib.pyplot as plt import numpy as np import seaborn seaborn.set_palette('husl') IMAGE_TEMPLATE = 'writeup_image%02d.png' def image0(): filename = IMAGE_TEMPLATE % (0,) num_points = 200 node_points = np.array([ [0, 5.0], [1.0, 9.0], [4.0, 14.0], ]) specific_s = np.array([ [1, 0, 0], [1, 0.5, 0.25], [1, 1, 1], ]) coeffs = np.linalg.solve(specific_s, node_points) s_values = np.linspace(0, 1, num_points) stacked_s = s_values[:, np.newaxis]**np.array([[0, 1, 2]]) x_y_vals = stacked_s.dot(coeffs) x_values = x_y_vals[:, 0] y_values = x_y_vals[:, 1] curved_edge, = plt.plot(x_values, y_values) edge_color = curved_edge.get_color() # Horizontal edge x_start = np.min(node_points[:, 0]) x_stop = np.max(node_points[:, 0]) y_start = np.min(node_points[:, 1]) plt.plot(np.linspace(x_start, x_stop, num_points), y_start * np.ones(num_points), color=edge_color) y_stop = np.max(node_points[:, 1]) # Vertical edge plt.plot(x_stop * np.ones(num_points), np.linspace(y_start, y_stop, num_points), color=edge_color) marked_x = node_points[:, 0] marked_y = node_points[:, 1] horizontal_midpt = 0.5 * (x_start + x_stop) vertical_midpt = 0.5 * (y_start + y_stop) marked_x = np.append(marked_x, [x_stop, horizontal_midpt, x_stop]) marked_y = np.append(marked_y, [y_start, y_start, vertical_midpt]) plt.plot(marked_x, marked_y, color='black', marker='o', linestyle='None') plt.axis('equal') x_width = x_stop - x_start y_width = y_stop - y_start plt.xlim(horizontal_midpt - 0.6 * x_width, horizontal_midpt + 0.6 * x_width) plt.ylim(vertical_midpt - 0.6 * y_width, vertical_midpt + 0.6 * y_width) plt.axis('off') plt.savefig(filename)
 import matplotlib.pyplot as plt import numpy as np import sympy def case1(x, y, s, f_symb): # CASE 1: ISOLATED POINT (-20, 12) AT s = (-1 +/- SQRT(-35)) / 6 # y1**2 = (6*x2 - 35)*(12*x2 + 35)**2 # x = 5*x2/2 - y1/3 - 305/24 # y = -3*x2/2 + y1/3 + 61/8 x_s = - 9 * s**3 + 18 * s**2 - 2 * s + 1 y_s = 9 * s**3 - 9 * s**2 + 5 * s f = (x**3 + 3 * x**2 * y + 3 * x * y**2 + y**3 - 5 * x**2 - 37 * x * y - 41 * y**2 + 52 * x + 52 * y - 48) # CHECK: f.subs({x: x_s, y: y_s}).expand() fx = f.diff(x) fy = f.diff(y) res_x = sympy.resultant(fx, fy, y) res_y = sympy.resultant(fx, fy, x) solns = sympy.solve([res_x, res_y]) for soln in solns: f_val = f.subs(soln) soln[f_symb] = f_val # f(-20, 12) = 0 # f_x(-20 , 12 ) = f_y(-20 , 12 ) = 0 # f_x(-20 , 13/4) = f_y(-20 , 13/4) = 0 # f_x(-65/12, 12 ) = f_y(-65/12, 12 ) = 0 # f_x(-65/12, 13/4) = f_y(-65/12, 13/4) = 0 # sympy.resultant(f, f.diff(y), y).factor() # sympy.solve(sympy.resultant(f, f.diff(y), y).factor()) # sympy.resultant(f, f.diff(x), x).factor() # sympy.solve(sympy.resultant(f, f.diff(x), x).factor()) return x_s, y_s, f, solns def case1_numerical(): x_verts = [-20, 5 - np.sqrt(4000.0/243), 5 + np.sqrt(4000.0/243)] y_verts = np.linspace(-10, 80, 500) ones_verts = np.ones(y_verts.shape) for x_vert in x_verts: curr_x = ones_verts * x_vert plt.plot(curr_x, y_verts, color='blue') x_for_y_vert = np.linspace(-40, 20, ones_verts.size) y_vert = 12 # No real roots in 243*y**2 - 486*y + 275 curr_y = y_vert * ones_verts plt.plot(x_for_y_vert, curr_y, color='red') min_x = -25 max_x = 10 num_pts = int((max_x - min_x) * 100) x_values = np.linspace(-25, 10, num_pts + 1) scatter_pts = [] for x in x_values: coeffs = [ 1, 3 * x - 41, 3 * x**2 - 37 * x + 52, x**3 - 5 * x**2 + 52 * x - 48, ] roots = np.roots(coeffs) roots = np.real_if_close(roots) # root, = roots[np.isreal(roots)] real_roots = roots[np.isreal(roots)] if len(real_roots) < 1: raise ValueError('Expected at least one real') # print x, len(real_roots) for root in real_roots: scatter_pts.append((x, root.real)) # root = real_roots[0] # root = root.real # y_values.append(root) # y_values = np.array(y_values) scatter_pts = np.array(scatter_pts) plt.scatter(scatter_pts[:, 0], scatter_pts[:, 1], s=3) plt.axis('equal') plt.show() def case2(x, y, s, f_symb): # CASE 2: SELF-INTERSECTION (-5399/1296, 41149/2592) AT s = (37 +/- SQRT(2229)) / 72 x_s = - 36 * s**3 + 36 * s**2 + 7 * s - 4 y_s = - 18 * s**3 + 72 * s**2 - 52 * s + 7 f = (2 * x**3 - 12 * x**2 * y + 24 * x * y**2 - 16 * y**3 + 1428 * x**2 + 660 * x * y + 744 * y**2 - 6320 * x - 5393 * y - 16689) / 2 # CHECK: f.subs({x: x_s, y: y_s}).expand() fx = f.diff(x) fy = f.diff(y) res_x = sympy.resultant(fx, fy, y) res_y = sympy.resultant(fx, fy, x) solns = sympy.solve([res_x, res_y]) for soln in solns: f_val = f.subs(soln) soln[f_symb] = f_val # FX = f.diff(x).subs({x: x_s, y: y_s}).factor() # FY = f.diff(y).subs({x: x_s, y: y_s}).factor() # EXTRA = (FY / x_s.diff(s)).factor() # CHECK: EXTRA.coeff(s**2) > 0 # CHECK: EXTRA + (FX / y_s.diff(s)).factor() == 0 return x_s, y_s, f, solns def case3(x, y, s, f_symb): # CASE 3: NO STATIONARY POINT (VALID CUBIC) x_s = s**3 + s**2 + s + 1 y_s = 2 * s**3 + 2 * s**2 + 1 f = (y - 2 * x)**3 + (y - 2 * x)**2 + 2 * x + 3 * y - 5 # CHECK: f.subs({x: x_s, y: y_s}).expand() fx = f.diff(x) fy = f.diff(y) res_x = sympy.resultant(fx, fy, y) res_y = sympy.resultant(fx, fy, x) solns = sympy.solve([res_x, res_y]) for soln in solns: f_val = f.subs(soln) soln[f_symb] = f_val return x_s, y_s, f, solns def case4(x, y, s, f_symb): # CASE 4: ISOLATED POINT (-20, -48) AT s = 1 +/- SQRT(-6) # f(-20 + h, -48 + k) # = 0 + 0 * h + 0 * k - 105 * h**2 + 96 * h * k - 22 * k**2 - # 8 * h**3 + 12 * h**2 * k - 6 * h * k**2 + k**3 # = kp**3 - 22 * kp**2 + 8 * h * kp - h**2 # WHERE kp = k - 2 * h x_s = s**3 + s**2 + s + 1 y_s = 2 * s**3 + 3 * s**2 + 1 f = (y - 2 * x)**3 - 9 * x**2 + 2 * y**2 + 24 * x - 16 # CHECK: f.subs({x: x_s, y: y_s}).expand() fx = f.diff(x) fy = f.diff(y) res_x = sympy.resultant(fx, fy, y) res_y = sympy.resultant(fx, fy, x) solns = sympy.solve([res_x, res_y]) for soln in solns: f_val = f.subs(soln) soln[f_symb] = f_val return x_s, y_s, f, solns def case5(x, y, s, f_symb): # CASE 5: ISOLATED POINT (307/125, 392/125) AT s = (-1 +/- SQRT(-3)) / 5 x_s = 4 * s**3 + 5 * s**2 + 2 * s + 3 y_s = - s**3 + 5 * s**2 + 2 * s + 4 f = ((x + 4 * y)**3 - 832 * x**2 + 344 * x * y - 937 * y**2 + 2333 * x + 2332 * y - 4834) # CHECK: f.subs({x: x_s, y: y_s}).expand() fx = f.diff(x) fy = f.diff(y) res_x = sympy.resultant(fx, fy, y) res_y = sympy.resultant(fx, fy, x) solns = sympy.solve([res_x, res_y]) for soln in solns: f_val = f.subs(soln) soln[f_symb] = f_val return x_s, y_s, f, solns def case6(x, y, s, f_symb): # CASE 6: SELF-INTERSECTION (-1/4, 11/8) AT s = (-1 +/- SQRT(5)) / 4 # s = 1/2 --> (x, y) = (-3/4, 13/8) # --> GRAD(f) = (-5, -14), dGAMMA/ds = (-7/2, 5/4) # --> DIR_DERIV(f, ( 5/sqrt(221), 14/sqrt(221) ) ) = -sqrt(221) # --> DIR_DERIV(f, (-5/sqrt(221), -14/sqrt(221) ) ) = sqrt(221) x_s = -2 * s**3 - 2 * s**2 y_s = - s**3 + s**2 + s + 1 f = ((x - 2 * y)**3 + 20 * x**2 - 24 * x * y + 32 * y**2 + 16 * x - 40 * y + 16) # CHECK: f.subs({x: x_s, y: y_s}).expand() fx = f.diff(x) fy = f.diff(y) res_x = sympy.resultant(fx, fy, y) res_y = sympy.resultant(fx, fy, x) solns = sympy.solve([res_x, res_y]) for soln in solns: f_val = f.subs(soln) soln[f_symb] = f_val return x_s, y_s, f, solns def case7(x, y, s, f_symb): # CASE 7: CUSP (0, 0) AT s = 0 x_s = s**2 y_s = s**3 f = y**2 - x**3 # CHECK: f.subs({x: x_s, y: y_s}).expand() fx = f.diff(x) fy = f.diff(y) res_x = sympy.resultant(fx, fy, y) res_y = sympy.resultant(fx, fy, x) solns = sympy.solve([res_x, res_y]) for soln in solns: f_val = f.subs(soln) soln[f_symb] = f_val return x_s, y_s, f, solns def case8(x, y, s, f_symb): # CASE 8: CUSP (-20, 8) AT s = 2 # x'(s) = 6(s + 1)(s - 2) # y'(s) = -(s - 2)(3s + 2) x_s = 2 * s**3 - 3 * s**2 - 12 * s y_s = - s**3 + 2 * s**2 + 4 * s f = (x + 2 * y)**3 - 4 * x**2 - 24 * x * y - 33 * y**2 - 16 * x - 48 * y # CHECK: f.subs({x: x_s, y: y_s}).expand() fx = f.diff(x) fy = f.diff(y) res_x = sympy.resultant(fx, fy, y) res_y = sympy.resultant(fx, fy, x) solns = sympy.solve([res_x, res_y]) for soln in solns: f_val = f.subs(soln) soln[f_symb] = f_val return x_s, y_s, f, solns def case8_plot(): s_vals = np.linspace(2 - 3, 2 + 1.5, 1000) x_vals = 2 * s_vals**3 - 3 * s_vals**2 - 12 * s_vals y_vals = - s_vals**3 + 2 * s_vals**2 + 4 * s_vals plt.plot(x_vals, y_vals) plt.show() def case9(x, y, s, f_symb): # CASE 9: CUSP (7, -4) AT s = -1 # x'(s) = 6(s + 1)(s - 2) # y'(s) = (s + 1)(7 - 3s) x_s = 2 * s**3 - 3 * s**2 - 12 * s y_s = - s**3 + 2 * s**2 + 7 * s f = (x + 2 * y)**3 - 22 * x**2 - 78 * x * y - 69 * y**2 - 7 * x - 12 * y # CHECK: f.subs({x: x_s, y: y_s}).expand() fx = f.diff(x) fy = f.diff(y) res_x = sympy.resultant(fx, fy, y) res_y = sympy.resultant(fx, fy, x) solns = sympy.solve([res_x, res_y]) for soln in solns: f_val = f.subs(soln) soln[f_symb] = f_val return x_s, y_s, f, solns def case9_plot(): s_vals = np.linspace(-1 - 2, -1 + 3.5, 1000) x_vals = 2 * s_vals**3 - 3 * s_vals**2 - 12 * s_vals y_vals = - s_vals**3 + 2 * s_vals**2 + 7 * s_vals plt.plot(x_vals, y_vals) plt.show() def main(): x, y, s, f_symb = sympy.symbols('x, y, s, f') x_s, y_s, f, solns = case1(x, y, s, f_symb) for soln in solns: print soln # g = 2 * s**2 + 1 # h = -s**3 + s**2 + 2 * s # SH = sympy.resultant(g - x, h - y, s) # sympy.resultant(g - x, g.diff(s), s) # sympy.resultant(h - y, h.diff(s), s) # sympy.Poly(SH, x) # sympy.Poly(SH, y) # SH = 8y^2 + (-x^3) + ... # sympy.solve((h - sympy.Rational(9, 8)).factor())