Chain of circles internally tangent to an ellipse. I tried to get an answer to this question (which was hastily closed) but couldn't find a proof, so I decided to ask it again, adding some of my efforts.
Suppose we have a finite sequence of $n$ circles ($n\ge10$, see figure below) whose centres lie on the major axis of an ellipse. All circles are internally tangent to the ellipse and each circle is also externally tangent to the preceding and following circle (if they exist). If $r_1$, $r_2$, ..., $r_n$ are the radii of these circles, prove that:
$$
r_7(r_1 + r_7) = r_4(r_4 + r_{10}).
$$

If $x_0$, $x_1$, ..., $x_n$ are the abscissae of the intersection points between the circles and the major axis (taking as origin the ellipse centre, see figure above), then it is not difficult to find a recursive relation for $x_k$. Let $a$, $b$ be the semi-major axis and semi-minor axis of the ellipse, $A$ and $B$ its foci, $O$ its centre and $c=AO=BO=\sqrt{a^2-b^2}$. If $C_k$ is the centre of $k$-th circle and $P_k$ one of its tangency points with the ellipse, then radius $P_kC_k$ is the normal to the ellipse at $P_k$ and thus the bisector of $\angle AP_kB$. It follows from the length of bisector formula that
$$
P_kC_k={b^2\over a^2}\sqrt{AP_k\cdot BP_k}={b^2\over c^2}\sqrt{AC_k\cdot BC_k}=
{b^2\over c^2}\sqrt{c^2-c_k^2},
$$
where $c_k$ is the abscissa of centre $C_k$. 
Inserting here $P_kC_k=(x_{k}-x_{k-1})/2$ and $c_k=(x_{k}+x_{k-1})/2$, 
then squaring both sides and rearranging, one finds:
$$
x_k^2+x_{k-1}^2-2(2e^2-1)x_kx_{k-1}=4e^2b^2,
$$
where $e=c/a$ is the eccentricity of the ellipse. From the above recursive equation one can find, once $x_0$ is given, all $x_k$ and thus compute $r_k=(x_{k}-x_{k-1})/2$ for all values of $k$. I used these   results with GeoGebra to draw the first figure, and could numerically check that the formula to prove holds for any value of $x_0$.
Nonetheless, I couldn't obtain a real proof of that formula using algebra, hence I believe I'm missing a simpler way to find those radii. Any idea to prove the statement is welcome.

 A: In the theory of Lucas sequences 
if $\,v\,$ is constant, and
$$ x_{n+1} = v\,x_n - x_{n-1} \tag{1}$$
 for all $\,n,\,(P=v,\;Q=1)\,$ then
$$ x_n^2 - x_{n+1}x_{n-1} = u \tag{2}$$
 for a constant $\,u.\,$ This implies that
$$ x_n^2 -v\,x_n x_{n+1} + x_{n+1}^2 = u \tag{3}$$
for all $\,n\,$ because
$$ x_{n-1}\!+\!x_{n+1} \!=\! v\,x_n, \;\;
x_{n-1}x_{n+1}\!=\!x_n^2\!-\!u. \tag{4}$$
In your case, the constants are
$$ u=4e^2b^2,\quad v=4e^2-2. \tag{5}$$
Also, check that if
$$ r_k:=(x_{k}-x_{k-1})/2, \tag{6} $$ then
$$ (r_n+r_{n+6})/r_{n+3} = v^3-3v \tag{7} $$
 for all $\,n.$ This
implies that
$$ r_{n+6}(r_n + r_{n+6}) = r_{n+3}(r_{n+3}+r_{n+6}). \tag{8}$$ for all $\,n.\,$ By the way, there is a
similar result for $\,(x_{n+m}+x_{n-m})/x_n\,$ and
$\,(r_{n+m}+r_{n-m})/r_n\,$ for any integer $\,m.$
Note that my answer is completely based on equation $(3)$
which was given in the question. I have not used any of
the geometric content of the question.
A: Let $P := (a \cos\phi, b \sin\phi)$ on an origin-centered ellipse with radii $a$ and $b$; define $c := \sqrt{a^2-b^2}$, so that the ellipse's eccentricity is $e := c/a$. The line through $P$, normal to the ellipse —that is, in direction $(b\cos\phi,a\sin\phi)$— meets the $x$-axis at $K:= (k,0)$, where $k:= c^2/a \cos\phi$. So, $K$ is the center of a circle internally tangent to the ellipse at $P$, and its radius, $r$, is given by
$$r^2 = |PK|^2 = \frac{b^2(a^2-c^2\cos^2\phi)}{a^2} = \frac{b^2(c^2-k^2)}{c^2} \tag{1}$$
so that
$$\frac{r^2}{b^2}+\frac{k^2}{c^2}=1 \tag{2}$$
This allow us to write, for some $\theta$,
$$r = b\sin\theta \qquad k = c \cos\theta \tag{3}$$
Now, suppose $\bigcirc K_0$ and $\bigcirc K_1$ are circles internally tangent to the ellipse, with respective centers and radii given by $(3)$ for $\theta = \theta_0$ and $\theta=\theta_1$. If these circles are tangent to each other (with $K_1$ "on the right" of $K_0$), then
$$\begin{align}
k_0 + r_0 &= k_1 - r_1 \\[4pt]
\to\quad -2 c \sin\frac{\theta_0 + \theta_1}{2} \sin\frac{\theta_0 - \theta_1}{2} &= -2 b \sin\frac{\theta_0 + \theta_1}{2} \cos\frac{\theta_0 - \theta_1}{2} \\[6pt]
\to\quad \tan\frac{\theta_0 - \theta_1}{2} &= \frac{b}{c} \\[6pt]
\to\quad \theta_1 &= \theta_0 - 2\arctan\frac{b}{c} \\[6pt]
&= \theta_0 - 2\arccos e  \tag{4}
\end{align}$$
More generally, if circles $\bigcirc K_i$, defined by $\theta = \theta_i$ in $(3)$, form a tangent chain, then
$$\theta_i = \theta_0 - 2 i \arccos e \tag{5}$$
where index $i$ is subject to certain viability conditions (eg, $\theta_i \geq 0$) that we'll assume hold. Thus, defining $\varepsilon := 2\arccos e$, we have
$$\begin{align}
\frac{r_{i+j} + r_{i-j}}{r_i} &= \frac{b\sin(\theta_0-(i+j)\psi)+b\sin(\theta_0-(i-j)\varepsilon)}{b \sin(\theta_0-i\varepsilon)} \\[6pt]
&= 2\cos j\varepsilon = 2\cos( 2j \arccos e ) \\[4pt]
&= 2\,T_{2j}(e) \tag{6}
\end{align}$$
where $T_{2j}$ is the $2j$-th Chebyshev polynomial of the first kind. Notably, the value of $(6)$ is independent of $i$. In particular, if we take $j=3$ and both $i=4$ and $i=7$, we can write
$$\frac{r_{4-3}+r_{4+3}}{r_4} = 2\;T_{2\cdot 3}(e) =\frac{r_{7-3}+r_{7+3}}{r_7} \tag{7}$$
which gives the result. $\square$

Addendum. In this follow-up question, @g.kov asks when an ellipse allows a "perfect packing" of $n$ tangent circles along its axis. It seems reasonable to append here a justification of the condition given there.
In a perfect packing, the first and last circles in a chain are tangent to the ellipse at the endpoints of the axis, so that their radii match the ellipse's radius of curvature (namely, $b^2/a$) at those points. Thus, we have
$$r_0 = r_{n-1} = \frac{b^2}{a} \quad\to\quad \sin\theta_0 = \sin\theta_{n-1} = \frac{b}{a} \quad\to\quad \cos\theta_0 = \cos\theta_{n-1} = e \tag{8}$$
We can say that $\theta_0 = \pi - \arccos e$ and $\theta_{n-1} = \arccos e$. By $(5)$, this implies
$$\arccos e = \theta_{n-1} = \theta_0 - 2(n-1)\arccos e = (\pi - \arccos e) - 2(n-1)\arccos e \tag{9}$$
so that
$$\pi = 2n\arccos e \qquad\to\qquad \cos \frac{\pi}{2n} = e \tag{10}$$
This is equivalent to @g.kov's condition for a perfectly-packable ellipse. $\square$
A: This isn't an answer, it's more of an annotation of the answer by @Somos, but it's too big to fit into the comment section.  The answer starts off by stating that if $\,v\,$ is constant, and
$$ x_{n+1} = v\,x_n - x_{n-1} \tag{1}$$
 for all $n,$ then
$$ x_n^2 - x_{n+1}x_{n-1} = u \tag{2}$$
 for a constant $\,u.\,$
(2) isn't obvious but can be derived from (1) by induction as follows.  Suppose (2) holds for a given value of $n$.  Then
$$
\begin{align}
x_{n+1}^2-x_{n+2}x_n & = x_{n+1}^2-(vx_{n+1}-x_n)x_n \\[1.2ex]
 & = x_{n+1}^2-vx_{n+1}x_n+x_n^2 \\[1.2ex] 
 & = x_{\mathstrut n}^2-x_{n+1}(vx_n-x_{n+1}) \\[1.2ex]
 & = x_{\mathstrut n}^2-x_{n+1}(vx_n-vx_n+x_{n-1}) \\[1.2ex] 
 & = x_{\mathstrut n}^2-x_{n+1}x_{n-1} \\[1.2ex]
 & = u.
\end{align}
$$
A: HINT:
In some sense the recursion of the radiuses is of order $1$.
Consider an ellipse of equation
$$x^2 + C^2 y^2 = 1$$ $C > 1$ parametrized as
$$(\cos t, \frac{1}{C} \sin t)$$
Now consider a point on the $x$ axis inside the ellipse, of coordinate $(x,0)$. If the point is not too close to the ends, then the closest distance to the ellipse is
$$r_x^2 = \frac{C^2(1-x^2)-1}{C^2(C^2-1)}$$
To see the limits for $x$, solve the equation
$$r_x^2 = (1-x)^2$$
(meaning that the closest point is in fact at the very end of the ellipse), with solution
$$x = 1- \frac{1}{C^2}$$
Therefore, the condition is $|x|\le 1-\frac{1}{C^2}(< 1)$
We certainly need a very flat ellipse ($C>>1$) for this problem to work.
