Scaling an Ellipse to be Tangent to Another Ellipse with a Different a/b Ratio

Question

I need to scale an ellipse to be tangent to another ellipse without moving them from their centers $(h_1,k_1)$ $(h_2,k_2)$ or modifying their a/b ratio.

$$\frac{(x-h)^2}{(a*s)^2} + \frac{(y-k)^2}{(b*s)^2} = 1$$

I'm looking for a fast analytical geometry solution to this problem, I understand that we could use numerical approach to solve this problem (e.g. Newtons Method)

If the a/b ratio of both ellipses were the same we could simply take advantage of the distance formula and radius at angle formula for ellipse to calculate a scale modifier and it would work in any situation and angle(demo of a working example):

$$ \frac{a \cdot b}{\sqrt{a^2 \sin^2 \theta + b^2 \cos^2 \theta}} $$

I believe this could be done by finding the two non-vertical common tangent points (NOT tangent lines) of the two ellipses. However, I couldn’t find a working solution anywhere. Since I’m using the solution in software, I need to be able to solve any possible equations with the quadratic formula.

In this demo, the scale of one of the ellipses is calculated by approximately finding the common tangent point on the first ellipse using the following formula:

$$ \begin{aligned} \theta &= \tan^{-1} \left( -\frac{b}{a} \cdot m \right) \\ x_1 &= h + a \cos \theta \\ y_1 &= k + b \sin \theta \\ x_2 &= h - a \cos \theta \\ y_2 &= k - b \sin \theta \\ \end{aligned} $$

Please let me know if you have any other questions or concerns.

Answer 1

If $r = \begin{bmatrix} x \\ y \end{bmatrix}$,

$C_1 = \begin{bmatrix} h_1 \\ k_1 \end{bmatrix} $

$ C_2 =\begin{bmatrix} h_2 \\ k_2 \end{bmatrix} $

$ Q_1 = \begin{bmatrix} \dfrac{1}{a_1^2} && 0 \\ 0 && \dfrac{1}{b_1^2} \end{bmatrix} $

$ Q_2 = \begin{bmatrix} \dfrac{1}{a_2^2} && 0 \\ 0 && \dfrac{1}{b_2^2} \end{bmatrix} $

Then the first ellipse is given by

$ (r - C_1)^T Q_1 (r - C_1) = 1 \tag{1} $

And the scaled second ellipse is given by

$ (r - C_2)^T Q_2 (r - C_2) = s^2 \tag{2}$

The gradient (normal vector) of the first ellipse is

$ \nabla_1 = 2 Q_1 (r - C_1 ) $

And the gradient of the second scaled ellipse is

$ \nabla_2 = 2 Q_2 (r - C_2 ) $

And we want

$ \nabla_1 = k \ \nabla_2 $

Define the matrix

$E = e_1 e_2^T - e_2 e_1^T = \begin{bmatrix} 0 && 1 \\ -1 && 0 \end{bmatrix} $

Then we require that

$ (r - C_1)^T Q_1^T E Q_2 (r - C_2) = 0 \tag{3}$

Equations $(1)$ and $(3)$ are a system of two quadratic equations in two unknowns, which are the $x$ and $y$ coordinates of the $r$, The tangency point that lies on the first (fixed) ellipse.

To solve them analytically, express the solutions of $(1)$ in parametric form, then you'll have

$ r = (x,y) = (h_1 + a_1 \cos t , k_1 + b_1 \sin t ) , \ t \in [0, 2 \pi) \tag{4}$

On the other hand, equation $(3)$ when expanded, we give

$ H x y + A x + B y + C = 0 \tag{5} $

where

$ H = \dfrac{1}{(a_1 b_2)^2} - \dfrac{1}{(a_2 b_1)^2} $

$ A = \dfrac{k_1}{(a_2 b_1)^2} - \dfrac{k_2}{(a_1 b_2)^2} $

$ B = \dfrac{h_2}{(a_2 b_1)^2} - \dfrac{h_1}{(a_1 b_2)^2} $

$ C = \dfrac{h_1 k_2}{(a_1 b_2)^2} - \dfrac{h_2 k_1}{(a_2 b_1)^2} $

Substituting $(4)$ into $(5)$ give us

$ H ( h_1 + a_1 \cos t )(k_1 + b_1 \sin t) + A (h_1 + a_1 \cos t) + B ( k_1 + b_1 \sin t ) + C = 0 $

Multiplying this we get the equation

$ A_1 \cos t + A_2 \sin t + A_3 \cos t \sin t + A_4 = 0 \tag{6}$

where

$ A_1 = A a_1 + H a_1 k_1 $

$ A_2 = B b_1 + H b_1 h_1 $

$ A_3 = H a_1 b_1 $

$ A_4 = H h_1 k_1 + A h_1 + B k_1 + C $

To solve equation $(6)$, we use the transformation

$ z = \tan \left( \dfrac{t}{2} \right) $

Then it would follow that

$ \cos t = \dfrac{1 - z^2}{1 + z^2} $ and $ \sin t =\dfrac{2 z }{1 + z^2} $

Plug these into $(6)$ and multiply through by $(1 + z^2)^2$. This gives you,

$ A_1 (1 - z^2) (1 + z^2) + A_2 (2 z)(1 + z^2) + A_3 (2 z)(1 - z^2) + A_4 (1 + z^2)^2 = 0 \tag{7}$

Expanding $(7)$ gives

$ a_4 z^4 + a_3 z^3 + a_2 z^2 + a_1 z + a_0 = 0 \tag{8} $

where

$a_4 = - A_1 + A_4 $

$a_3 = 2 A_2 - 2 A_3 $

$a_2 = 2 A_4 $

$a_1 = 2 A_2 + 2 A_3 $

$a_0 = A_1 + A_4 $

Feeding this quartic polynomial (quartic means of degree 4) to a library routine that finds polynomial roots generates all the $4$ roots, out of which, we're only interested in the real roots. There should be exactly two real roots. You can write your own quartic polynomial solver. See here for details. I've included my routine that I use to find the roots of a quartic polynomial below.

Once a root $z$ are found, the corresponding $t$ is calculated using $t = 2 \tan^{-1} z$, and from there, we can calculate $r$ from equation $(4)$. Once $r$ is found, the value of $s$ can be calculated from equation $(2)$.

The above analysis applies only to the case where the ellipses have their axes parallel to the coordinate axes. But the quadratic equations $(1), (2), (3)$ apply to any two ellipses in general orientation. In this case, we can still apply the same procedure, but the equations are far more complex. However, the two equations $(1)$ and $(3)$ can be solved directly for $r$ by a general two-variable quadratic equation solver.

Numerical Example:

Let $C_1 = (10, 0) $, $C_2 = (0, 6) $, $a_1 = 3 , b_1 = 2 $, $a_2 = 2, b_2 = 1 $. Feeding the above quadratic system to a solver generates the following solutions. The two values of $s$ are

$s_1 = 6.01335710947491 $

and

$ s_2 = 9.63595852170716 $

$s_1$ corresponds to external tangency, and $s_2$ corresponds to internal tangency.

Below are the two solutions (external and internal tangency)

The system of three quadratic equations was obtained for a pair of rotated ellipses. The first ellipse, has its center at $C_1 = (8,0)$, with $a_1 = 5, b_1 = 3$, and rotation angle of $30^\circ$. The second ellipse is centered at $C_2 = (0, 4)$, with $a_2 = 4 , b_2 = 2 $ and a rotation angle of $60^\circ$.

The scale factors obtained were

$s_1 = 2.679569 $ (for external tangency)

and

$ s_2 = 6.3 $ (for internal tangency)

They are plotted below.

---

---

Public Sub quartic1(a1, b1, c1, d1, e1, r)
Dim rt(4, 2) As Double
Dim a, b, c, d, e As Double
Dim i, j, k As Integer
Dim cubic_roots(3) As Double
Dim r0 As Double

' a z^4 + b z^3 + c z^2 + d z + e = 0

' normalizing


a = a1
b = b1
c = c1
d = d1
e = e1


If Abs(a) < 0.00000001 Then
  MsgBox ("This is not a quartic polynomial.  Invalid Entry") 
  Exit Sub
End If


b = b / a
c = c / a
d = d / a
e = e / a
a = 1


'  z^4 + b z^3 + c z^2 + d z + e = 0

' solve the resolvent cubic

a2 = -c
a1 = b * d - 4 * e
a0 = 4 * c * e - d ^ 2 - b ^ 2 * e

Call find_cubic_real_root(a2, a1, a0, x0)

r0 = b ^ 2 / 4 - c + x0

If Abs(r0) < 0.00000001 Then r0 = 0

If r0 >= 0 Then

   r0x = Sqr(r0)
   r0y = 0

Else

   r0x = 0
   r0y = Sqr(-r0)

End If

r0norm = Sqr(r0x ^ 2 + r0y ^ 2)

r2x = r0x ^ 2 - r0y ^ 2
r2y = 0

If Abs(r0norm) > 0.00000001 Then
   rinvx = 1 / r0norm ^ 2 * (r0x)
   rinvy = 1 / r0norm ^ 2 * (-r0y)

   d2x = 0.75 * b ^ 2 - r2x - 2 * c + 0.25 * (4 * b * c - 8 * d - b ^ 3) * rinvx
   d2y = 0.25 * (4 * b * c - 8 * d - b ^ 3) * rinvy

 ' now find the square root of d2

  theta = atn2(d2x, d2y)

   d2norm = Sqr(d2x ^ 2 + d2y ^ 2)

   d3x = d2norm ^ (1 / 2) * Cos(theta / 2)
   d3y = d2norm ^ (1 / 2) * Sin(theta / 2)

Else
   x1 = x0 ^ 2 - 4 * e

   If Abs(x1) < 0.00000001 Then x1 = 0
     d2 = 0.75 * b ^ 2 - 2 * c + 2 * Sqr(x1)
     If d2 >= 0 Then
        d3x = Sqr(d2)
        d3y = 0

     Else
        d3x = 0
        d3y = Sqr(-d2)
     End If

 End If


rt(1, 1) = -b / 4 + 0.5 * r0x + 0.5 * d3x
rt(1, 2) = 0.5 * r0y + 0.5 * d3y

rt(2, 1) = -b / 4 + 0.5 * r0x - 0.5 * d3x
rt(2, 2) = 0.5 * r0y - 0.5 * d3y


lnext_couple:

If r0norm > 0.00000001 Then
rinvx = 1 / r0norm ^ 2 * (r0x)
rinvy = 1 / r0norm ^ 2 * (-r0y)

e2x = 0.75 * b ^ 2 - r2x - 2 * c - 0.25 * (4 * b * c - 8 * d - b ^ 3) * rinvx
e2y = -0.25 * (4 * b * c - 8 * d - b ^ 3) * rinvy


theta = atn2(e2x, e2y)

e2norm = Sqr(e2x ^ 2 + e2y ^ 2)

e3x = e2norm ^ (1 / 2) * Cos(theta / 2)
e3y = e2norm ^ (1 / 2) * Sin(theta / 2)


Else
x1 = x0 ^ 2 - 4 * e
If Abs(x1) < 0.00000001 Then x1 = 0
e2 = 0.75 * b ^ 2 - 2 * c - 2 * Sqr(x1)

If e2 >= 0 Then

e3x = Sqr(e2)
e3y = 0

Else

e3x = 0
e3y = Sqr(-e2)

End If

End If

rt(3, 1) = -b / 4 - 0.5 * r0x + 0.5 * e3x
rt(3, 2) = -0.5 * r0y + 0.5 * e3y

rt(4, 1) = -b / 4 - 0.5 * r0x - 0.5 * e3x
rt(4, 2) = -0.5 * r0y - 0.5 * e3y

iroot = 0

For i = 1 To 4

If Abs(rt(i, 2)) < 0.000001 Then

  iroot = iroot + 1
  r(iroot) = rt(i, 1)

End If


lnext_i:
Next i

r(0) = iroot

' order roots in ascending order

For i = 1 To r(0) - 1
x0min = r(i)
jmin = i
For j = i + 1 To r(0)

If r(j) < x0min Then
x0min = r(j)
jmin = j
End If

Next j

'exchange i,j
If i <> jmin Then
   t = r(i)
   r(i) = r(jmin)
   r(jmin) = t
End If

Next i

End Sub

---

Public Sub find_cubic_real_root(a2, a1, a0, x0)
Dim b As Double
Dim q0 As Double

q = (3 * a1 - a2 ^ 2) / 9
r = (9 * a2 * a1 - 27 * a0 - 2 * a2 ^ 3) / 54

d = q ^ 3 + r ^ 2

p = WorksheetFunction.pi()

If d >= 0 Then
   q0 = r + Sqr(d)
   s = Sgn(q0) * (Abs(r + Sqr(d))) ^ (1 / 3)
   q0 = r - Sqr(d)
   t = Sgn(q0) * (Abs(r - Sqr(d))) ^ (1 / 3)
   x0 = -a2 / 3 + s + t

Else

   dsi = Sqr(-d)

   theta = WorkSheetFunction.Atan2(r, dsi)

   x0 = -a2 / 3 + 2 * (Sqr(r ^ 2 + dsi ^ 2)) ^ (1 / 3) * Cos(theta / 3)

End If


End Sub

---

Answer 2

Given two tangent ellipses

$$ \cases{ \frac{(x-h_1)^2}{a^2}+\frac{(y-k_1)^2}{b^2}=s^2\\ \frac{(x-h_2)^2}{a^2}+\frac{(y-k_2)^2}{b^2}=1 } $$

substituting $x$ we have

$$ 4a(b^2(h_1-h_2)^2+a^2(k_1-k_2)^2)y^2+-4 a( b^2 (h_1 - h_2)^2 (k_1 + k_2) + a^2 (k_1 - k_2) (b^2 + k_1^2 - k_2^2 - b^2 s^2))y+(b^4 (h1 - h2)^4)/a + a^3 (k_2^2 - k_1^2 + b^2 (s^2-1))^2 - 2 a b^2 (h_1 - h_2)^2 (b^2 (1 + s^2)-k_1^2 - k_2^2)=0 $$

or

$$ c_2y^2+c_1y+c_0=0 $$

but by tangency proprieties we know that

$$ c_2y^2+c_1y+c_0=d_0(y-y_0)^2,\ \ \forall y $$

and as a consequence

$$ \cases{ c_0-d_0 y_0^2=0\\ c_1+2d_0 y_0 =0\\ c_2 - d_0 = 0 } $$

so we have three equations and three unknowns $(s,d_0,y_0)$. Solving for $s$ we obtain

$$ s = \pm 1\pm \frac{1}{ab}\sqrt{a^2(k_1-k_2)^2+b^2(h_1-h_2)^2} $$

four solutions. Two considering $s^2$

Answer 3

This is what I had in mind; the attractive feature is that it results in a polynomial of degree four in one auxiliary variable I'm calling $t.$ The relevant value for $t$ may be found by your own numerical methods. The time I remember I had to numerically solve a one-variable thing, I id it by bisection. I recommend you do this yourself, not some canned solver.

What I did was move one ellipse to center around the origin, then scaling in just one direction is used, to make the other one a circle, centered at some $(d,e).$ The ellipse is now some (I use $c^2$ here) $$ \frac{x^2}{a^2} + \frac{y^2}{b^2} = c^2 $$ The gradient of the ellipse constraint function is $$ \left( \frac{2x}{a^2} , \frac{2y}{b^2} \right) $$ We have an unknown point $(x,y)$ on the ellipse and know the normal vector there. We need the point plus some real $t$ times the normal to equal the fixed center $(d,e)$ or $$ (x,y)+ t \left( \frac{2x}{a^2} , \frac{2y}{b^2} \right) = (d,e) $$ Well, $$ x \left( 1 + \frac{2t}{a^2} \right) = d $$ $$ y \left( 1 + \frac{2t}{b^2} \right) = e $$

We can solve in terms of $t$ $$ x = \frac{d}{ \left( 1 + \frac{2t}{a^2} \right) } = \frac{a^2d}{a^2 + 2t} $$

$$ y = \frac{e}{ \left( 1 + \frac{2t}{b^2} \right) } = \frac{b^2e}{b^2 + 2t} $$

Next $$ \frac{x^2}{a^2} = \frac{a^2d^2}{(a^2 + 2t)^2} $$ while $$ \frac{y^2}{b^2} = \frac{b^2e^2}{(b^2 + 2t)^2} $$

So $$ c^2 = \frac{a^2d^2}{(a^2 + 2t)^2} + \frac{b^2e^2}{(b^2 + 2t)^2} $$

with $$ c^2 (2t+a^2)^2 (2t+b^2)^2 = a^2 d^2 (2t+b^2)^2 + b^2 e^2 (2t+a^2)^2 $$

which does have a $t^4$ term on the left hand side as well as $t^3$ terms. I would start constructing my solver with numerous error reports along the way, essentially asking for aspects of a graph of the quartic. One you know exactly what root you need, continue with bisection or finish with Newton's. But there are occasional surprises with Newton's