# Finding Equations for Common Tangents Between Two Ellipses, A General Solution

I'm interested in finding the four common tangents between two ellipses. While I've found some fascinating approaches using dual conics to identify their intersection points (link 1, link 2), my primary objective is to obtain the tangent equations in the original Cartesian coordinate system.

The methods I've come across seem to work within the realm of projective geometry, which is beyond my current mathematical expertise. Any help in simplifying this to Cartesian coordinates would be invaluable.

#### My Current Understanding

Let's say we have two ellipses defined by their canonical matrices E1 and E2.

1. Dual Conics: The common tangents can be found at the intersection points of the dual conics C1 = E1^-1 and C2 = E2^-1.

2. Intersection Points: To find these, we consider the pencil of conics C1 - t * C2 (where t is real) and solve for t in det(C1 - t * C2) = 0. Since this is a cubic equation, it yields up to three roots t0, t1, and t2, which correspond to up to three degenerate conics—usually intersecting lines. The equation for these degenerate conics is C1 - t_i * C2 = 0.

3. Coordinates of Intersection Points: As these points lie on the conics, for each t_i we should have:

(C1 - t_i * C2) * X_i = 0


Subsequently:

((C2^-1 * C1) - t_i * I) * X_i = 0


Based on an answer by @Futurologist here, this is an eigenvalue equation, with t_i as the eigenvalue and X_i as the eigenvector. Therefore, X_i should give me the coordinates of the intersection points, which should also provide the coefficients for the tangents in Cartesian coordinates.

Have I understood this correctly, or is there something I'm missing?

Example Let's have a concrete example with the 2 ellipsis:

E1: 36 * x^2 + 9 * y^2 - 324 = 0

and

E2: 4 * x^2 + y^2 - 56 * x - 6 * y + 201 = 0

Considering the method I described above, please can you explain how you get the four tangents final equations?

• How are your ellipses given? Commented Oct 17, 2023 at 13:54
• @trula I've added more details to my question
– MohG
Commented Oct 17, 2023 at 15:14
• Maybe not what you're looking for but here's a treatment by an old master: Salmon, Conic Sections, pg 344 Article 378 Commented Oct 17, 2023 at 20:59
• Note that there are possible from 4 to none such tangents.
– Moti
Commented Oct 18, 2023 at 4:56
• A real affine ellipse is the zero-set of $\frac{x^2}{a^2}+\frac{y^2}{b^2}-1.$ All ellipses are generated from these ones by rotation and translation or is the zero-set of $$\text{Constant}\cdot(\frac{(\cos{\theta_i}(x-h_i)+\sin{\theta_i}(y-k_i))^2}{a^2}+\frac{(-\sin{\theta_i}(x-h_i)+\cos{\theta_i}(y-k_i))^2}{b^2}-1), i=1,2$$ Commented Oct 18, 2023 at 9:53

Update: complex numbers

This implementation with complex numbers seems to work for now

import math
import cmath
import numpy as np

'''
Algorithm for finding the coefficients of
the linear equations of the four tangent
lines of a pair of non-intersecting ellipses
'''

# given a pair of conics, as a pair of symmetric matrices,
# calculate the vector k = (k[0], k[1], k[2]) of values for each of which
# the conic c1 - k[i]*c2 from the pencil of conics c1 - t*c2
# is a degenerate conic (the anti-symmetric product of a pair of linear forms)
# and also find the matrix U
# of the projective transformation that simplifies the geometry of
# the pair of conics, and the geometry of the pencil c1 - t*c2 in general,
# as well as the geometry of the three degenerate conics in particular

def get_transformation(c1, c2):
'''
c1 and c2 are 3 by 3 symmetric matrices of the two conics
'''
c21 = np.linalg.inv(c2).dot(c1)
k, U = np.linalg.eig(c21)
return k, U

def find_common_points(c1, c2):
k, U = get_transformation(c1, c2)
L1 = (U.T).dot((c1 - k[0]*c2).dot(U))
L2 = (U.T).dot((c1 - k[1]*c2).dot(U))
x_0 = cmath.sqrt(-(L2[2,2] / L2[0,0]))
y_0 = cmath.sqrt(-(L1[2,2] / L1[1,1]))
sol = np.array([[ x_0,  y_0, 1],
[-x_0, -y_0, 1],
[-x_0,  y_0, 1],
[ x_0, -y_0, 1]])
sol = sol.dot(U.T)
return sol

# combination of the algorithms form above with the
# transformation of the problem of common tangents of the two ellipses
# to its dual - the intersection points of the two dual ellipses
# The intersection points of the dual ellipses
# are the common tangents of the two original ellipses

def find_common_tangents(c1, c2):
dc1 = np.linalg.inv(c1)
dc2 = np.linalg.inv(c2)
return find_common_points(dc1, dc2)


Previous version: real version, non-intersecting ellipses

If an algorithm is what you are seeking, here is one, based on my previous posts you cited:

import math
import numpy as np

'''
Algorithm for finding the coefficients of
the linear equations of the four tangent
lines of a pair of non-intersecting ellipses
'''

# given a pair of conics, as a pair of symmetric matrices,
# calculate the vector k = (k[0], k[1], k[2]) of values for each of which
# the conic c1 - k[i]*c2 from the pencil of conics c1 - t*c2
# is a degenerate conic (the anti-symmetric product of a pair of linear forms)
# and also find the matrix U
# of the projective transformation that simplifies the geometry of
# the pair of conics, and the geometry of the pencil c1 - t*c2 in general,
# as well as the geometry of the three degenerate conics in particular

def get_transformation(c1, c2):
'''
c1 and c2 are 3 by 3 symmetric matrices of the two conics
'''
c21 = np.linalg.inv(c2).dot(c1)
k, U = np.linalg.eig(c21)
return k, U

# find the common points, i.e. points of intersection,
# of a pair of non-intersecting ellipses
# represented by a pair of symmetric matrices

def find_common_points(c1, c2):
'''
c1 and c2 are 3 by 3 symmetric matrices of the two conics
'''
k, U = get_transformation(c1, c2)
L1 = (U.T).dot((c1 - k[0]*c2).dot(U))
L2 = (U.T).dot((c1 - k[1]*c2).dot(U))
sol = np.ones((4,3), dtype=float)
for i in range(2):
for j in range(2):
sol[i+2*j,0:2] = np.array([math.sqrt(abs(L2[2,2] / L2[0,0]))*(-1)**i, math.sqrt(abs(L1[2,2] / L1[1,1]))*(-1)**j])
sol = sol.dot(U.T)
return sol

# combination of the algorithms form above with the
# transformation of the problem of common tangents of the two ellipses
# to its dual - the intersection points of the two dual ellipses
# The intersection points of the dual ellipses
# are the common tangents of the two original ellipses

def find_common_tangents(c1, c2):
dc1 = np.linalg.inv(c1)
dc2 = np.linalg.inv(c2)
return find_common_points(dc1, dc2)


I tested it with one example, and it worked, but have not run more tests. Here is a test script:

'''
Begin generation of test conics
'''

def cos_sin(angle_deg):
return math.cos(angle_deg*math.pi/180), math.sin(angle_deg*math.pi/180)

def rotation(cs_sn):
return np.array([[cs_sn[0], -cs_sn[1]],
[cs_sn[1],  cs_sn[0]]])

def isom(angle, translation):
'''
#isometry from conic-aligned coordinate system (conic attached)
#to global coordinate system (world system)
isometry from global coordinate system (world system)
to conic-aligned coordinate system (conic attached)
'''
cos_, sin_ = cos_sin(-angle)
tr = - rotation((cos_, sin_)).dot(translation)
return np.array([[ cos_, -sin_, tr[0]],
[ sin_,  cos_, tr[1]],
[    0,     0,    1 ]])

def conic(major, minor, angle, center):
D = np.array([[minor**2,        0,                 0],
[       0, major**2,                 0],
[       0,        0, -(minor*major)**2]])
U = isom(angle, center)
return (U.T).dot(D.dot(U))

def conic_coeff(c):
return np.array( (c[0,0], 2*c[0,1], c[1,1], 2*c[0,2], 2*c[1,2], c[2,2]) )

'''End generation of test conics'''

# BEGIN test:

rnd=5 # rounding to 5 digits after the decimal point. Choose as you like.

a = 2
b = 1
cntr = np.array([0,0])
w = 45

C1 = conic(a, b, w, cntr)
Eq1 = conic_coeff(C1).round(rnd)

a = 3
b = 1
cntr = np.array([7,0])
w = 120

C2 = conic(a, b, w, cntr)
Eq2 = conic_coeff(C2).round(rnd)

common_tangents = find_common_tangents(C1, C2)

print()
print('equation of the conic 1:')
print(f'{Eq1[0]}x**2 + {Eq1[1]}xy + {Eq1[2]}y**2 + {Eq1[3]}x + {Eq1[4]}y + {Eq1[5]} = 0')
print()
print('equation of conic 2:')
print(f'{Eq2[0]}x**2 + {Eq2[1]}xy + {Eq2[2]}y**2 + {Eq2[3]}x + {Eq2[4]}y + {Eq2[5]} = 0')
print()

print('The coefficients of the linear equations of the four tangent lines to the conics 1 and 2 are: ')
print(common_tangents.round(rnd))


After running this particular test, I got as output:

equation of conic 1:
2.5x**2 + -3.0xy + 2.5y**2 + 0.0x + 0.0y + -4.0 = 0

equation of conic 2:
7.0x**2 + 6.9282xy + 3.0y**2 + -98.0x + -48.49742y + 334.0 = 0

the coefficients of the linear equations of the four tangent lines to the conics 1 and 2 are:
[[ 0.52106  0.85509 -1.96045]
[ 0.07325  0.63866  1.08326]
[-0.274    1.24082 -1.73691]
[-0.72181  1.02439  1.3068 ]]


Then I went to Desmos, because I find it convenient, and plotted the output, just to check visually:

• Amazing! You're a boss! One extra question: is there a simple way to pick only external tangents?
– MohG
Commented Oct 19, 2023 at 9:05
• And a second question: what about the cases when the 2 initial ellipses are intersecting. The script doesn't seem to find the right tangents
– MohG
Commented Oct 19, 2023 at 13:43
• @MohG 1) There is a way to find the inner common tangents but it requires some work: e.g. the centres of the ellipses should be on opposite sides of a common tangent. But your post started with: "[I am] interested in finding the four common tangents between two ellipses" 2) If the two ellipses intersect, then you have to modify the algorithm to work with complex numbers. Commented Oct 19, 2023 at 14:03
• @MohG In your OP you state: "Therefore, X_i should give me the coordinates of the intersection points, which should also provide the coefficients for the tangents in Cartesian coordinates." This is not true. The three eigenvectors X_i are the columns of the orthogonal matrix U, which is the projective transformation that transforms the two ellipses from general configuration into a special configuration in which the solutions to the problem are very easy to find. Then, one uses U to transform back the solutions of the special configuration into the solutions of the general configuration. Commented Oct 19, 2023 at 14:21
• @MohG It doesn't work for intersecting ellipses, because intersecting ellipses do not have inner tangents. Only outer tangents. If the ellipses intersect in two (real) points, there are only two real tangents (the other two are complex) and they are outer. When the ellipses intersect in four (real) points, then there are four real tangents and they are all outer. Commented Nov 15, 2023 at 16:05

You propose E0: $$324\cdot(x^2/3^2+y^2/6^2-1) =0$$ and E1: $$4\cdot((x-7)^2+(y-3)^2/2^2-1)=0$$ which have duals, dropping the constants $$324$$ and $$4$$; $$3^2X^2+6^2Y^2-1=0$$ (the ones centered on the origin are particularly easy) and $$12X^2+\frac{21}2 XY+\frac54 Y^2+\frac72 X+\frac32 Y+\frac14=0.$$ The matrix of the second ellipse is $$M=\begin{pmatrix}1&0&-7\\0&\frac14&-\frac34\\-7&-\frac34&\frac{201}4\end{pmatrix}$$ and that has inverse $$\frac1{\det M}\begin{pmatrix}12&\frac{21}4&\frac74\\\frac{21}4&\frac54&\frac34\\\frac74&\frac34&\frac14\end{pmatrix}.$$ The ideal of the two duals has primary decomposition $$\langle 21X+9Y+2,205Y^2+4Y-5\rangle \cap \langle 21X+9Y+4,615Y^2+24Y-11\rangle$$ so they intersect in $$(X,Y)=$$ $$(-\frac{9\sqrt3\sqrt7+56}{615},\frac{7\sqrt3\sqrt7-2}{205}), (\frac{9\sqrt3\sqrt7-56}{615}, -\frac{7\sqrt3\sqrt7+2}{205}),\\(-\frac{3\sqrt3\sqrt{47}+112}{615}, \frac{7 \sqrt3 \sqrt{47}-12}{615}),(\frac{3 \sqrt3\sqrt{47}-112}{615}, -\frac{7\sqrt3\sqrt{47}+12}{615}),$$ or $$(-0.1581189939099228,0.1467220969009311),\\(-0.02399482722828852,-0.1662342920228823),\\ (-0.2400374410749817,0.1156429180638462), \\(-0.124190201201441,-0.1546673083077486),$$ see geogebra, making the four common tangents $$-\frac{9\sqrt3\sqrt7+56}{615}x+\frac{7\sqrt3\sqrt7-2}{205}y+1=0, \frac{9\sqrt3\sqrt7-56}{615}x -\frac{7\sqrt3\sqrt7+2}{205}y+1=0\\-\frac{3\sqrt3\sqrt{47}+112}{615}x +\frac{7 \sqrt3 \sqrt{47}-12}{615}y+1=0,\frac{3 \sqrt3\sqrt{47}-112}{615}x -\frac{7\sqrt3\sqrt{47}+12}{615}y+1=0$$ or $$-0.1581189939099228x+0.1467220969009311y+1=0,\\-0.02399482722828852x-0.1662342920228823y+1=0,\\ -0.2400374410749817x+0.1156429180638462y+1=0, \\-0.124190201201441x-0.1546673083077486y+1=0.$$

The affine geometry is in the projective one you've seen in other posts as $$z=1$$ and dually $$Z=1,$$ the universal line then goes from being $$Xx+Yx+Zz=0$$ to $$Xx+Yx+1=0.$$

Edit

If one ellipse is in standard form $$\frac{x^2}{a_1^2}+\frac{y^2}{b_1^2}-1=0,$$ which you can achieve by rotating and translating the plane, landing one in standard form and the other in general form $$\frac{(\cos{\alpha}(x-h)+\sin{\alpha}(y-k))^2}{a_2^2}+\frac{(-\sin{\alpha}(x-h)+\cos{\alpha}(y-k))^2}{b_2^2}-1=0.$$ Solving this system, I proceed by producing a Lex grobner basis: the first two elements of which are given below. The first is a quartic in $$y:$$ $$q_4(y)=0$$ and the next is of the form $$Cx+p_3(y)=0.$$ Solving the quartic you get one $$x$$-value for each real root of $$q_4(y)=0.$$ But there's no way in general around solving a quartic equation for this problem.

(a1^4*a2^4*sin(α)^4-2*a1^2*a2^4*b1^2*sin(α)^4+a2^4*b1^4*sin(α)^4-2*a1^4*a2^2*b2^2*sin(α)^4+4*a1^2*a2^2*b1^2*b2^2*sin(α)^4-2*a2^2*b1^4*b2^2*sin(α)^4+a1^4*b2^4*sin(α)^4-2*a1^2*b1^2*b2^4*sin(α)^4+b1^4*b2^4*sin(α)^4-8*a1^2*a2^2*b1^2*h*k*cos(α)*sin(α)+8*a1^2*b1^2*b2^2*h*k*cos(α)*sin(α)+2*a1^2*a2^4*b1^2*sin(α)^2-2*a2^4*b1^4*sin(α)^2+2*a1^4*a2^2*b2^2*sin(α)^2-4*a1^2*a2^2*b1^2*b2^2*sin(α)^2+2*a2^2*b1^4*b2^2*sin(α)^2-2*a1^4*b2^4*sin(α)^2+2*a1^2*b1^2*b2^4*sin(α)^2+2*a1^2*a2^2*b1^2*h^2*sin(α)^2+2*a2^2*b1^4*h^2*sin(α)^2-2*a1^2*b1^2*b2^2*h^2*sin(α)^2-2*b1^4*b2^2*h^2*sin(α)^2-2*a1^4*a2^2*k^2*sin(α)^2-2*a1^2*a2^2*b1^2*k^2*sin(α)^2+2*a1^4*b2^2*k^2*sin(α)^2+2*a1^2*b1^2*b2^2*k^2*sin(α)^2+a2^4*b1^4-2*a1^2*a2^2*b1^2*b2^2+a1^4*b2^4-2*a2^2*b1^4*h^2+2*a1^2*b1^2*b2^2*h^2+b1^4*h^4+2*a1^2*a2^2*b1^2*k^2-2*a1^4*b2^2*k^2+2*a1^2*b1^2*h^2*k^2+a1^4*k^4)*y^4+(-8*a1^2*a2^2*b1^2*h*cos(α)*sin(α)+8*a1^2*b1^2*b2^2*h*cos(α)*sin(α)-4*a1^4*a2^2*k*sin(α)^2-4*a1^2*a2^2*b1^2*k*sin(α)^2+4*a1^4*b2^2*k*sin(α)^2+4*a1^2*b1^2*b2^2*k*sin(α)^2+4*a1^2*a2^2*b1^2*k-4*a1^4*b2^2*k+4*a1^2*b1^2*h^2*k+4*a1^4*k^3)*y^3+(2*a1^2*a2^4*sin(α)^4-2*a2^4*b1^2*sin(α)^4-4*a1^2*a2^2*b2^2*sin(α)^4+4*a2^2*b1^2*b2^2*sin(α)^4+2*a1^2*b2^4*sin(α)^4-2*b1^2*b2^4*sin(α)^4+8*a1^2*a2^2*h*k*cos(α)*sin(α)-8*a1^2*b2^2*h*k*cos(α)*sin(α)-2*a1^4*a2^2*sin(α)^2-2*a1^2*a2^4*sin(α)^2-2*a1^2*a2^2*b1^2*sin(α)^2+4*a2^4*b1^2*sin(α)^2+2*a1^4*b2^2*sin(α)^2+4*a1^2*a2^2*b2^2*sin(α)^2+2*a1^2*b1^2*b2^2*sin(α)^2-4*a2^2*b1^2*b2^2*sin(α)^2-2*a1^2*b2^4*sin(α)^2-2*a1^2*a2^2*h^2*sin(α)^2-4*a2^2*b1^2*h^2*sin(α)^2+2*a1^2*b2^2*h^2*sin(α)^2+4*b1^2*b2^2*h^2*sin(α)^2+2*a1^2*a2^2*k^2*sin(α)^2-2*a1^2*b2^2*k^2*sin(α)^2+2*a1^2*a2^2*b1^2-2*a2^4*b1^2-2*a1^4*b2^2+2*a1^2*a2^2*b2^2+2*a1^2*b1^2*h^2+4*a2^2*b1^2*h^2-2*a1^2*b2^2*h^2-2*b1^2*h^4+6*a1^4*k^2-2*a1^2*a2^2*k^2-2*a1^2*h^2*k^2)*y^2+(8*a1^2*a2^2*h*cos(α)*sin(α)-8*a1^2*b2^2*h*cos(α)*sin(α)+4*a1^2*a2^2*k*sin(α)^2-4*a1^2*b2^2*k*sin(α)^2+4*a1^4*k-4*a1^2*a2^2*k-4*a1^2*h^2*k)*y+a2^4*sin(α)^4-2*a2^2*b2^2*sin(α)^4+b2^4*sin(α)^4+2*a1^2*a2^2*sin(α)^2-2*a2^4*sin(α)^2-2*a1^2*b2^2*sin(α)^2+2*a2^2*b2^2*sin(α)^2+2*a2^2*h^2*sin(α)^2-2*b2^2*h^2*sin(α)^2+a1^4-2*a1^2*a2^2+a2^4-2*a1^2*h^2-2*a2^2*h^2+h^4=0

(2*a1^2*a2^6*sin(α)^6-6*a1^2*a2^4*b2^2*sin(α)^6+6*a1^2*a2^2*b2^4*sin(α)^6-2*a1^2*b2^6*sin(α)^6+4*a1^2*a2^4*h*k*cos(α)*sin(α)^3-8*a1^2*a2^2*b2^2*h*k*cos(α)*sin(α)^3+4*a1^2*b2^4*h*k*cos(α)*sin(α)^3+2*a1^4*a2^4*sin(α)^4-4*a1^2*a2^6*sin(α)^4-4*a1^4*a2^2*b2^2*sin(α)^4+10*a1^2*a2^4*b2^2*sin(α)^4+2*a1^4*b2^4*sin(α)^4-8*a1^2*a2^2*b2^4*sin(α)^4+2*a1^2*b2^6*sin(α)^4+2*a1^2*a2^4*h^2*sin(α)^4-4*a1^2*a2^2*b2^2*h^2*sin(α)^4+2*a1^2*b2^4*h^2*sin(α)^4-4*a1^2*a2^4*h*k*cos(α)*sin(α)+4*a1^2*a2^2*b2^2*h*k*cos(α)*sin(α)+4*a1^2*a2^2*h^3*k*cos(α)*sin(α)-4*a1^2*b2^2*h^3*k*cos(α)*sin(α)-2*a1^4*a2^4*sin(α)^2+2*a1^2*a2^6*sin(α)^2+4*a1^4*a2^2*b2^2*sin(α)^2-4*a1^2*a2^4*b2^2*sin(α)^2-2*a1^4*b2^4*sin(α)^2+2*a1^2*a2^2*b2^4*sin(α)^2+2*a1^4*a2^2*h^2*sin(α)^2-2*a1^2*a2^4*h^2*sin(α)^2+2*a1^2*a2^2*b1^2*h^2*sin(α)^2-2*a1^4*b2^2*h^2*sin(α)^2+4*a1^2*a2^2*b2^2*h^2*sin(α)^2-2*a1^2*b1^2*b2^2*h^2*sin(α)^2-2*a1^2*b2^4*h^2*sin(α)^2-2*a1^2*a2^2*h^2*k^2*sin(α)^2+2*a1^2*b2^2*h^2*k^2*sin(α)^2-2*a1^2*a2^2*b1^2*h^2+2*a1^4*b2^2*h^2+2*a1^2*b1^2*h^4+2*a1^2*a2^2*h^2*k^2-2*a1^2*h^4*k^2)*x+(-a1^4*a2^6*cos(α)*sin(α)^5+2*a1^2*a2^6*b1^2*cos(α)*sin(α)^5-a2^6*b1^4*cos(α)*sin(α)^5+3*a1^4*a2^4*b2^2*cos(α)*sin(α)^5-6*a1^2*a2^4*b1^2*b2^2*cos(α)*sin(α)^5+3*a2^4*b1^4*b2^2*cos(α)*sin(α)^5-3*a1^4*a2^2*b2^4*cos(α)*sin(α)^5+6*a1^2*a2^2*b1^2*b2^4*cos(α)*sin(α)^5-3*a2^2*b1^4*b2^4*cos(α)*sin(α)^5+a1^4*b2^6*cos(α)*sin(α)^5-2*a1^2*b1^2*b2^6*cos(α)*sin(α)^5+b1^4*b2^6*cos(α)*sin(α)^5-2*a1^2*a2^6*b1^2*cos(α)*sin(α)^3+2*a2^6*b1^4*cos(α)*sin(α)^3-2*a1^4*a2^4*b2^2*cos(α)*sin(α)^3+6*a1^2*a2^4*b1^2*b2^2*cos(α)*sin(α)^3-4*a2^4*b1^4*b2^2*cos(α)*sin(α)^3+4*a1^4*a2^2*b2^4*cos(α)*sin(α)^3-6*a1^2*a2^2*b1^2*b2^4*cos(α)*sin(α)^3+2*a2^2*b1^4*b2^4*cos(α)*sin(α)^3-2*a1^4*b2^6*cos(α)*sin(α)^3+2*a1^2*b1^2*b2^6*cos(α)*sin(α)^3-2*a1^2*a2^4*b1^2*h^2*cos(α)*sin(α)^3-2*a2^4*b1^4*h^2*cos(α)*sin(α)^3+4*a1^2*a2^2*b1^2*b2^2*h^2*cos(α)*sin(α)^3+4*a2^2*b1^4*b2^2*h^2*cos(α)*sin(α)^3-2*a1^2*b1^2*b2^4*h^2*cos(α)*sin(α)^3-2*b1^4*b2^4*h^2*cos(α)*sin(α)^3+2*a1^4*a2^4*k^2*cos(α)*sin(α)^3+2*a1^2*a2^4*b1^2*k^2*cos(α)*sin(α)^3-4*a1^4*a2^2*b2^2*k^2*cos(α)*sin(α)^3-4*a1^2*a2^2*b1^2*b2^2*k^2*cos(α)*sin(α)^3+2*a1^4*b2^4*k^2*cos(α)*sin(α)^3+2*a1^2*b1^2*b2^4*k^2*cos(α)*sin(α)^3+a1^4*a2^4*h*k*sin(α)^4-10*a1^2*a2^4*b1^2*h*k*sin(α)^4+a2^4*b1^4*h*k*sin(α)^4-2*a1^4*a2^2*b2^2*h*k*sin(α)^4+20*a1^2*a2^2*b1^2*b2^2*h*k*sin(α)^4-2*a2^2*b1^4*b2^2*h*k*sin(α)^4+a1^4*b2^4*h*k*sin(α)^4-10*a1^2*b1^2*b2^4*h*k*sin(α)^4+b1^4*b2^4*h*k*sin(α)^4-a2^6*b1^4*cos(α)*sin(α)+2*a1^2*a2^4*b1^2*b2^2*cos(α)*sin(α)+a2^4*b1^4*b2^2*cos(α)*sin(α)-a1^4*a2^2*b2^4*cos(α)*sin(α)-2*a1^2*a2^2*b1^2*b2^4*cos(α)*sin(α)+a1^4*b2^6*cos(α)*sin(α)+2*a2^4*b1^4*h^2*cos(α)*sin(α)-2*a1^2*a2^2*b1^2*b2^2*h^2*cos(α)*sin(α)-2*a2^2*b1^4*b2^2*h^2*cos(α)*sin(α)+2*a1^2*b1^2*b2^4*h^2*cos(α)*sin(α)-a2^2*b1^4*h^4*cos(α)*sin(α)+b1^4*b2^2*h^4*cos(α)*sin(α)-2*a1^2*a2^4*b1^2*k^2*cos(α)*sin(α)+2*a1^4*a2^2*b2^2*k^2*cos(α)*sin(α)+2*a1^2*a2^2*b1^2*b2^2*k^2*cos(α)*sin(α)-2*a1^4*b2^4*k^2*cos(α)*sin(α)-10*a1^2*a2^2*b1^2*h^2*k^2*cos(α)*sin(α)+10*a1^2*b1^2*b2^2*h^2*k^2*cos(α)*sin(α)-a1^4*a2^2*k^4*cos(α)*sin(α)+a1^4*b2^2*k^4*cos(α)*sin(α)+10*a1^2*a2^4*b1^2*h*k*sin(α)^2-2*a2^4*b1^4*h*k*sin(α)^2+2*a1^4*a2^2*b2^2*h*k*sin(α)^2-20*a1^2*a2^2*b1^2*b2^2*h*k*sin(α)^2+2*a2^2*b1^4*b2^2*h*k*sin(α)^2-2*a1^4*b2^4*h*k*sin(α)^2+10*a1^2*b1^2*b2^4*h*k*sin(α)^2+2*a1^2*a2^2*b1^2*h^3*k*sin(α)^2+2*a2^2*b1^4*h^3*k*sin(α)^2-2*a1^2*b1^2*b2^2*h^3*k*sin(α)^2-2*b1^4*b2^2*h^3*k*sin(α)^2-2*a1^4*a2^2*h*k^3*sin(α)^2-2*a1^2*a2^2*b1^2*h*k^3*sin(α)^2+2*a1^4*b2^2*h*k^3*sin(α)^2+2*a1^2*b1^2*b2^2*h*k^3*sin(α)^2+a2^4*b1^4*h*k-2*a1^2*a2^2*b1^2*b2^2*h*k+a1^4*b2^4*h*k-2*a2^2*b1^4*h^3*k+2*a1^2*b1^2*b2^2*h^3*k+b1^4*h^5*k+2*a1^2*a2^2*b1^2*h*k^3-2*a1^4*b2^2*h*k^3+2*a1^2*b1^2*h^3*k^3+a1^4*h*k^5)*y^3+(4*a1^4*a2^4*k*cos(α)*sin(α)^3+4*a1^2*a2^4*b1^2*k*cos(α)*sin(α)^3-8*a1^4*a2^2*b2^2*k*cos(α)*sin(α)^3-8*a1^2*a2^2*b1^2*b2^2*k*cos(α)*sin(α)^3+4*a1^4*b2^4*k*cos(α)*sin(α)^3+4*a1^2*b1^2*b2^4*k*cos(α)*sin(α)^3-a1^4*a2^4*h*sin(α)^4-6*a1^2*a2^4*b1^2*h*sin(α)^4-a2^4*b1^4*h*sin(α)^4+2*a1^4*a2^2*b2^2*h*sin(α)^4+12*a1^2*a2^2*b1^2*b2^2*h*sin(α)^4+2*a2^2*b1^4*b2^2*h*sin(α)^4-a1^4*b2^4*h*sin(α)^4-6*a1^2*b1^2*b2^4*h*sin(α)^4-b1^4*b2^4*h*sin(α)^4-4*a1^2*a2^4*b1^2*k*cos(α)*sin(α)+4*a1^4*a2^2*b2^2*k*cos(α)*sin(α)+4*a1^2*a2^2*b1^2*b2^2*k*cos(α)*sin(α)-4*a1^4*b2^4*k*cos(α)*sin(α)-4*a1^2*a2^2*b1^2*h^2*k*cos(α)*sin(α)+4*a1^2*b1^2*b2^2*h^2*k*cos(α)*sin(α)-4*a1^4*a2^2*k^3*cos(α)*sin(α)+4*a1^4*b2^2*k^3*cos(α)*sin(α)+6*a1^2*a2^4*b1^2*h*sin(α)^2+2*a2^4*b1^4*h*sin(α)^2-2*a1^4*a2^2*b2^2*h*sin(α)^2-12*a1^2*a2^2*b1^2*b2^2*h*sin(α)^2-2*a2^2*b1^4*b2^2*h*sin(α)^2+2*a1^4*b2^4*h*sin(α)^2+6*a1^2*b1^2*b2^4*h*sin(α)^2-2*a1^2*a2^2*b1^2*h^3*sin(α)^2-2*a2^2*b1^4*h^3*sin(α)^2+2*a1^2*b1^2*b2^2*h^3*sin(α)^2+2*b1^4*b2^2*h^3*sin(α)^2-2*a1^4*a2^2*h*k^2*sin(α)^2-2*a1^2*a2^2*b1^2*h*k^2*sin(α)^2+2*a1^4*b2^2*h*k^2*sin(α)^2+2*a1^2*b1^2*b2^2*h*k^2*sin(α)^2-a2^4*b1^4*h+2*a1^2*a2^2*b1^2*b2^2*h-a1^4*b2^4*h+2*a2^2*b1^4*h^3-2*a1^2*b1^2*b2^2*h^3-b1^4*h^5+2*a1^2*a2^2*b1^2*h*k^2-2*a1^4*b2^2*h*k^2+2*a1^2*b1^2*h^3*k^2+3*a1^4*h*k^4)*y^2+(-3*a1^2*a2^6*cos(α)*sin(α)^5+a2^6*b1^2*cos(α)*sin(α)^5+9*a1^2*a2^4*b2^2*cos(α)*sin(α)^5-3*a2^4*b1^2*b2^2*cos(α)*sin(α)^5-9*a1^2*a2^2*b2^4*cos(α)*sin(α)^5+3*a2^2*b1^2*b2^4*cos(α)*sin(α)^5+3*a1^2*b2^6*cos(α)*sin(α)^5-b1^2*b2^6*cos(α)*sin(α)^5+a1^4*a2^4*cos(α)*sin(α)^3+3*a1^2*a2^6*cos(α)*sin(α)^3+a1^2*a2^4*b1^2*cos(α)*sin(α)^3-2*a2^6*b1^2*cos(α)*sin(α)^3-2*a1^4*a2^2*b2^2*cos(α)*sin(α)^3-9*a1^2*a2^4*b2^2*cos(α)*sin(α)^3-2*a1^2*a2^2*b1^2*b2^2*cos(α)*sin(α)^3+4*a2^4*b1^2*b2^2*cos(α)*sin(α)^3+a1^4*b2^4*cos(α)*sin(α)^3+9*a1^2*a2^2*b2^4*cos(α)*sin(α)^3+a1^2*b1^2*b2^4*cos(α)*sin(α)^3-2*a2^2*b1^2*b2^4*cos(α)*sin(α)^3-3*a1^2*b2^6*cos(α)*sin(α)^3+a1^2*a2^4*h^2*cos(α)*sin(α)^3+2*a2^4*b1^2*h^2*cos(α)*sin(α)^3-2*a1^2*a2^2*b2^2*h^2*cos(α)*sin(α)^3-4*a2^2*b1^2*b2^2*h^2*cos(α)*sin(α)^3+a1^2*b2^4*h^2*cos(α)*sin(α)^3+2*b1^2*b2^4*h^2*cos(α)*sin(α)^3-a1^2*a2^4*k^2*cos(α)*sin(α)^3+2*a1^2*a2^2*b2^2*k^2*cos(α)*sin(α)^3-a1^2*b2^4*k^2*cos(α)*sin(α)^3+11*a1^2*a2^4*h*k*sin(α)^4-a2^4*b1^2*h*k*sin(α)^4-22*a1^2*a2^2*b2^2*h*k*sin(α)^4+2*a2^2*b1^2*b2^2*h*k*sin(α)^4+11*a1^2*b2^4*h*k*sin(α)^4-b1^2*b2^4*h*k*sin(α)^4-a1^2*a2^4*b1^2*cos(α)*sin(α)+a2^6*b1^2*cos(α)*sin(α)+a1^4*a2^2*b2^2*cos(α)*sin(α)-a1^2*a2^4*b2^2*cos(α)*sin(α)+a1^2*a2^2*b1^2*b2^2*cos(α)*sin(α)-a2^4*b1^2*b2^2*cos(α)*sin(α)-a1^4*b2^4*cos(α)*sin(α)+a1^2*a2^2*b2^4*cos(α)*sin(α)+a1^2*a2^2*b1^2*h^2*cos(α)*sin(α)-2*a2^4*b1^2*h^2*cos(α)*sin(α)+a1^2*a2^2*b2^2*h^2*cos(α)*sin(α)-a1^2*b1^2*b2^2*h^2*cos(α)*sin(α)+2*a2^2*b1^2*b2^2*h^2*cos(α)*sin(α)-a1^2*b2^4*h^2*cos(α)*sin(α)+a2^2*b1^2*h^4*cos(α)*sin(α)-b1^2*b2^2*h^4*cos(α)*sin(α)-5*a1^4*a2^2*k^2*cos(α)*sin(α)+a1^2*a2^4*k^2*cos(α)*sin(α)+5*a1^4*b2^2*k^2*cos(α)*sin(α)-a1^2*a2^2*b2^2*k^2*cos(α)*sin(α)+11*a1^2*a2^2*h^2*k^2*cos(α)*sin(α)-11*a1^2*b2^2*h^2*k^2*cos(α)*sin(α)+a1^4*a2^2*h*k*sin(α)^2-11*a1^2*a2^4*h*k*sin(α)^2+a1^2*a2^2*b1^2*h*k*sin(α)^2+2*a2^4*b1^2*h*k*sin(α)^2-a1^4*b2^2*h*k*sin(α)^2+22*a1^2*a2^2*b2^2*h*k*sin(α)^2-a1^2*b1^2*b2^2*h*k*sin(α)^2-2*a2^2*b1^2*b2^2*h*k*sin(α)^2-11*a1^2*b2^4*h*k*sin(α)^2-a1^2*a2^2*h^3*k*sin(α)^2-2*a2^2*b1^2*h^3*k*sin(α)^2+a1^2*b2^2*h^3*k*sin(α)^2+2*b1^2*b2^2*h^3*k*sin(α)^2+a1^2*a2^2*h*k^3*sin(α)^2-a1^2*b2^2*h*k^3*sin(α)^2-a1^2*a2^2*b1^2*h*k-a2^4*b1^2*h*k+a1^4*b2^2*h*k+a1^2*a2^2*b2^2*h*k+a1^2*b1^2*h^3*k+2*a2^2*b1^2*h^3*k-a1^2*b2^2*h^3*k-b1^2*h^5*k+3*a1^4*h*k^3-a1^2*a2^2*h*k^3-3*a1^2*h^3*k^3)*y-2*a1^2*a2^4*k*cos(α)*sin(α)^3+4*a1^2*a2^2*b2^2*k*cos(α)*sin(α)^3-2*a1^2*b2^4*k*cos(α)*sin(α)^3+5*a1^2*a2^4*h*sin(α)^4+a2^4*b1^2*h*sin(α)^4-10*a1^2*a2^2*b2^2*h*sin(α)^4-2*a2^2*b1^2*b2^2*h*sin(α)^4+5*a1^2*b2^4*h*sin(α)^4+b1^2*b2^4*h*sin(α)^4-2*a1^4*a2^2*k*cos(α)*sin(α)+2*a1^2*a2^4*k*cos(α)*sin(α)+2*a1^4*b2^2*k*cos(α)*sin(α)-2*a1^2*a2^2*b2^2*k*cos(α)*sin(α)+6*a1^2*a2^2*h^2*k*cos(α)*sin(α)-6*a1^2*b2^2*h^2*k*cos(α)*sin(α)+a1^4*a2^2*h*sin(α)^2-5*a1^2*a2^4*h*sin(α)^2+a1^2*a2^2*b1^2*h*sin(α)^2-2*a2^4*b1^2*h*sin(α)^2-a1^4*b2^2*h*sin(α)^2+10*a1^2*a2^2*b2^2*h*sin(α)^2-a1^2*b1^2*b2^2*h*sin(α)^2+2*a2^2*b1^2*b2^2*h*sin(α)^2-5*a1^2*b2^4*h*sin(α)^2+a1^2*a2^2*h^3*sin(α)^2+2*a2^2*b1^2*h^3*sin(α)^2-a1^2*b2^2*h^3*sin(α)^2-2*b1^2*b2^2*h^3*sin(α)^2+a1^2*a2^2*h*k^2*sin(α)^2-a1^2*b2^2*h*k^2*sin(α)^2-a1^2*a2^2*b1^2*h+a2^4*b1^2*h+a1^4*b2^2*h-a1^2*a2^2*b2^2*h+a1^2*b1^2*h^3-2*a2^2*b1^2*h^3+a1^2*b2^2*h^3+b1^2*h^5+a1^4*h*k^2-a1^2*a2^2*h*k^2-3*a1^2*h^3*k^2=0

• amazing but not sure to understand what is M and how do you compute the intersection points of the duals?
– MohG
Commented Oct 18, 2023 at 18:01
• Ok I think I understand how do you get the intersection points: you "just" solved manually the equation system of the 2 dual conics. But let say this equation is more complexe (with an x*y element for instance). Actually I'm looking for a method to get the intersections points "automatically" from the 2 duals conics equations. Something that I could implement in my software.
– MohG
Commented Oct 18, 2023 at 19:24
• @MohG Can you farm out the calculation to a CAS? Then that’s easy. The general solution to the form in my comment (with five parameters for each ellipse) to the question ran for a looong time without producing a solution like the one you seek. Commented Oct 18, 2023 at 19:42
• Nope. I would like to find a general solution without any expensive numerical computation. A solution has been posted here [stackoverflow.com/questions/65260321/… in python but I can’t make it work because not fully understood. Did you check?
– MohG
Commented Oct 18, 2023 at 20:10
• @MohG CASs are symbolic, which is better than numeric for conics. Also the CAS solution I alluded to in the previous comment would be easy to program (substitute the parameters) alas it takes way long to run with ten of these. If you have more information, it could lighten the computation, though. Commented Oct 18, 2023 at 20:17

Considering the ellipses

$$\cases{ 36x^2 + 9y^2 - 324=0\\ 4x^2 + y^2 - 56x - 6y + 201=0 }$$

and a generic line given by $$y = mx + n$$ if the ellipses are tangent to the line, then after substitution, the resulting polynomials in $$x$$ verify for all $$x$$:

$$\cases{ 36x^2 + 9(mx+n)^2 - 324=c_1(x-x_1)^2\\ 4x^2 + (mx +n)^2 - 56x - 6(mx+n) + 201=c_2(x-x_2)^2 }$$

which gives us the equations

$$\left\{ \begin{array} 9 n^2-c_1 x_1^2-324 =0\\ 2 c_1 x_1+18 m n =0\\ 9 m^2+36 -c_1=0\\ n^2-6 n+201 -c_2 x_2^2=0 \\ 2 c_2 x_2+2 m n-6 m-56=0 \\ m^2+4-c_2=0 \\ \end{array} \right.$$

Eliminating $$\{c_1,c_2,x_1,x_2\}$$ we arrive at

$$\cases{ 126075 - 7380 n - 5234 n^2 + 164 n^3 + 55 n^4=0\\ m = \frac{1}{6216}(14349 + 1435 n - 329 n^2 - 55 n^3) }$$

and after solving we get

$$\left[ \begin{array}{cc} m & n \\ \frac{1}{15} \left(7+2 \sqrt{21}\right) & -\frac{1}{5} \left(2+7 \sqrt{21}\right) \\ \frac{1}{15} \left(7-2 \sqrt{21}\right) & \frac{1}{5} \left(7 \sqrt{21}-2\right) \\ \frac{1}{33} \left(21+4 \sqrt{141}\right) & -\frac{1}{11} \left(12+7 \sqrt{141}\right) \\ \frac{1}{33} \left(21-4 \sqrt{141}\right) & \frac{1}{11} \left(7 \sqrt{141}-12\right) \\ \end{array} \right]$$