share/contrib/conics_04.mac

   1 /* Intersection of two conics
   2     David Billinghurst, 2021-09-02
   3   */
   4
   5 display2d:false;
   6
   7 "intersection of two arbitrary conics"$
   8
   9 eq1: a1*x^2+b1*x*y+c1*y^2+d1*x+e1*y+f1$
  10 eq2: a2*x^2+b2*x*y+c2*y^2+d2*x+e2*y+f2$
  11
  12 "Eliminate y^2"$
  13
  14 eq3a:c2*eq1-c1*eq2$
  15 eq3a:expand(eq3a);
  16
  17 "New variables to simplify."$
  18 subs:[
  19   a3=a1*c2-a2*c1,
  20   b3=b1*c2-b2*c1,
  21   d3=d1*c2-d2*c1,
  22   e3=e1*c2-e2*c1,
  23   f3=f1*c2-f2*c1];
  24
  25 eq3: a3*x^2+b3*x*y+d3*x+e3*y+f3;
  26
  27 "Show (eq3a) is equal to (eq3)"$
  28 expand(eq3-eq3a),subs;
  29
  30 "Solve eq3 for y
  31
  32 This step can't be performed when b_3*x + e_3 = 0.  This is not uncommon
  33 For example, it will happen if both conics are centred on the x-axis.
  34
  35 If b3=0 and e3=0 then (eq3) is already free of y and can be solved for x.
  36 When it occurs we should be able to substitute our solution for x into (eq1)
  37 and (eq2) and solve for y.  There may be spurious solutions that do not
  38 satisfy both equations.
  39
  40 The case x=0 and e_3=0 should be investigated."$
  41
  42 eq4:solve(eq3,y);
  43 eq4:eq4[1];
  44
  45 "Substitute (eq4) into (eq1).  Multiply through by denominator"$
  46
  47 eq5:subst(eq4,eq1);
  48 eq5*(b3*x+e3)^2;
  49 eq6:ratsimp(%);
  50
  51 "eq6 is a quartic in x, k4*x^4+k3*x^3+k2*x^2+k1*x+k0"$
  52 k4:coeff(eq6,x,4);
  53 k3:coeff(eq6,x,3);
  54 k2:coeff(eq6,x,2);
  55 k1:coeff(eq6,x,1);
  56 k0:coeff(eq6,x,0);
  57 "Confirm eq6-k4*x^4+k3*x^3+k2*x^2+k1*x+k0 = 0"$
  58 expand(eq6-(k4*x^4+k3*x^3+k2*x^2+k1*x+k0));
  59
  60 "Example 1"$
  61
  62 "Ellipse centred at origin with semi-major axis 2 and 1"$
  63 ellipse:(x/2)^2+y^2-1;
  64 ellipse:expand(ellipse);
  65 "Create a list of coefficients for later use"$
  66 coeffs_1:[a1=1/4,b1=0,c1=1,d1=0,e1=0,f1=-1];
  67
  68 "Circle centre (3,0) radius 2"$
  69 circle:(x-3)^2+y^2-2^2;
  70 circle:expand(circle);
  71 coeffs_2:[a2=1,b2=0,c2=1,d2=-6,e2=0,f2=5];
  72
  73 "Generate a list of coefficients to substitute"$
  74 coeffs:subs,coeffs_1,coeffs_2;
  75 coeffs:append(coeffs,coeffs_1,coeffs_2);
  76 s0:subst(coeffs,k0);
  77 s1:subst(coeffs,k1);
  78 s2:subst(coeffs,k2);
  79 s3:subst(coeffs,k3);
  80 s4:subst(coeffs,k4);
  81 p:s4*x^4+s3*x^3+s2*x^2+s1*x+s0;
  82
  83 "Solve quartic p and select real roots.  Only one (multiple) root"$
  84 soln:allroots(p);
  85 x1:rhs(soln[1]);
  86
  87 "Now solve for y.
  88
  89 Could try eq4: y = -(a3*x^2+d3*x+f3)/(b3*x+e3), but in this case b3=0
  90 and e3=0 so divide by zero
  91
  92 Alternatively substitute x into either conic eqn and solve for y, then
  93 check the
  94 solutions as some may be spurious."$
  95
  96 "First real root x1"$
  97 subst(x=x1,ellipse);
  98 ysolns:allroots(%);
  99 y1:rhs(ysolns[1]);
 100 y2:rhs(ysolns[2]);
 101 "Check if roots satisfy other equation circle"$
 102 circle,x=x1,y=y1;
 103 circle,x=x1,y=y2;
 104
 105 kill(ellipse,x1,y1,y2,coeffs_1,coeffs,p,s0,s1,s2,s3,s4);
 106
 107
 108 "Example 2:
 109
 110   Ellipse centred at origin
 111     semi-major axis length 2 in direction x=-y
 112     semi-minor axis length 1 in direction x=y
 113
 114   Same circle centre (3,0) radius 2"$
 115
 116 ellipse:(5*y^2)/8+(3*x*y)/4+(5*x^2)/8-1;
 117 ellipse:expand(ellipse);
 118 coeffs_1:[a1=5/8,b1=3/4,c1=5/8,d1=0,e1=0,f1=-1];
 119
 120 coeffs:subs,coeffs_1,coeffs_2;
 121 coeffs:append(coeffs,coeffs_1,coeffs_2);
 122 s0:k0,coeffs;
 123 s1:k1,coeffs;
 124 s2:k2,coeffs;
 125 s3:k3,coeffs;
 126 s4:k4,coeffs;
 127 p:s4*x^4+s3*x^3+s2*x^2+s1*x+s0;
 128
 129 soln:allroots(p);
 130 "Select the real roots"$
 131 x1:rhs(soln[1]);
 132 x2:rhs(soln[2]);
 133
 134 "Use eq4 as b3*x+e3 /= 0"$
 135 y1:rhs(eq4),x=x1,coeffs;
 136 y2:rhs(eq4),x=x2,coeffs;
 137 ellipse,x=x1,y=y1;
 138 circle,x=x1,y=y1;
 139 ellipse,x=x2,y=y2;
 140 circle,x=x2,y=y2;
 141
 142 "Or substitute x into conic equations, solve for y and select common roots"$
 143 allroots(subst(x=x1,circle));
 144 allroots(subst(x=x1,ellipse));
 145 allroots(subst(x=x2,circle));
 146 allroots(subst(x=x2,ellipse));