Cubic polynomial smoothly connecting two circles

Question

Given two circles with radii $R_L$ and $R_R$ and centers at $(-(R_L+a),\,0)$ and $(R_R+a,\,0)$, respectively, find a cubic polynomial $p(x)=b+cx^2+dx^3$ that smoothly connects the two circles.

$b$ is a parameter so $p(0)=b$ and the linear term of the polynomial is omitted because we want $\frac{\mathrm{d}p}{\mathrm{d}x}\Big|_{x=0}=0$.

My attempt at a solution. Let $\mathrm{C}_{L,R}$ the equations for the upper half of the $L,R$ circles. I formulate two equations relating $\mathrm{C}_{L,R}$ and $p$, and two equations relating $\frac{\mathrm{d}}{\mathrm{d}x}\mathrm{C}_{L,R}$ and $\frac{\mathrm{d}p}{\mathrm{d}x}$. Let $x_{L,R}$ be the points where $p(x)$ and $\mathrm{C}_{L,R}(x)$ intersect, then:

$$ \mathrm{C}_L(x_L)-p(x_L) = 0 $$ $$ \mathrm{C}_R(x_R)-p(x_R) = 0 $$ $$ \frac{\mathrm{d}}{\mathrm{d}x}\mathrm{C}_L(x_L) - \frac{\mathrm{d}p}{\mathrm{d}x}(x_L) = 0 $$ $$ \frac{\mathrm{d}}{\mathrm{d}x}\mathrm{C}_R(x_R) - \frac{\mathrm{d}p}{\mathrm{d}x}(x_R) = 0 $$ so I have a system of four nonlinear equations with four unknowns $(x_L,\,x_R,\,c,\,d)$. I coded a simple Newton's method for the system and it works well for some combinations of parameters $(a,\,b,\,R_L,\,R_R)$ when the initial guess is close enough, especially when $|R_L-R_R|$ is not too large and I use a constant damping for the Newton iterations. I can find initial guesses that I think are good via a graphical interface I coded. However, as $|R_L-R_R|$ gets larger the solver fails spectacularly to converge even with very close initial guesses and very small damping. (I should add that I'm actually taking the square of the equations to avoid square roots of negative numbers during the Newton iterations).

My question is threefold:

a) what other method or modification can I use to make the solver more stable?

b) this problem seems to me like it should be solved somewhere, do you know a reference?

c) more generally, is there a reason this should fail as horribly as it does when $|R_L-R_R|>>1$?

Answer 1

Here is a method giving all solutions in a deterministic way.

Take a look at the following figure :

Fig. 1 : The exhaustive set of 8 solution curves in the case $a=1, \ b=1, \ R_L=3, \ R_R=1$. Please note that some tangencies are internal to the circles.

How has this result been obtained ? By using a Computer Algebra System with a rather simple system of 4 polynomial equations in the 4 unknowns $c,d,x_L,x_R$ where (your notations) the two last parameters are the abscissas of the tangency points of the cubic curve with equation

$$y=f(x)=b+cx^2+dx^3$$

with the Left and Right circle resp.

Here they are:

$$\begin{cases}(x_L+R_L+a)^2+f(x_L)^2&=&R_L^2& \ \ (i)\\ (x_R-R_R-a)^2+f(x_R)^2&=&R_R^2& \ \ (ii)\\ \dfrac{f(x_L)}{x_L+R_L+a}&=&- \dfrac{1}{f'(x_L)}& \ \ (iii)\\ \dfrac{f(x_R)}{x_R-R_R-a}&=&- \dfrac{1}{f'(x_R)}& \ \ (iv)\\ \end{cases}\tag{1}$$

The first two equations express that $M_L=(x_L,f(x_L))$ and $M_R=(x_R,f(x_R))$ belong to their resp. circles.

Let us now explain the third equation. Let $C_L(-(R_L+a),0)$ denotes the center of the first circle. The slope of the tangent in $M_L$ is $f'(x_L)$. $\vec{C_LM_L}=\binom{x_L+R_L+a}{f(x_L)-0}$, being orthogonal to this tangent, has a slope $-\dfrac{1}{f'(x_L)}$ (see there). Similar reasoning for the fourth equation.

Remarks :

1) I haven't found the limitations you mention in the case of a large gap between $R_L$ and $R_R$. For example when $R_L=100$ and $R_R=1$, (with $a=1$ and $b=0$), one finds 14 solutions...

2) The third equation can be written : $f(x_L)f'(x_L)+x_L+R_L+a=0$ : this avoids possible division by zero. Same thing for the fourth equation.

Another case (Fig. 2), this time with $R_L=R_R$ displaying some spurious solutions under the form of ... second degree curves, i.e., parabolas (one can indeed have $d=0$) :

Fig. 2.

Still another case, with four solutions (one of them looks to be a parabola but isn't):

Fig. 3.

Here is the Matlab program that has given figure 1
(running time : around 30 seconds on my rather slow computer). Please note the "isreal" conditions (we want the different unknowns to be real) : in fact in the first case, there are $40$ solutions, $32$ of them being spurious, i.e., with complex coefficients...

Final remark : In fact, there exists a different way to solve system (1) by doing it in two steps ; first, by expressing the fact that equations (i) and (iii) have a common root $x_L$ giving (using a "resultant" ) a first (non-linear) constraint between $c$ and $d$. Doing the same for equations (ii) and (iv) for common root $x_R$, gives a second constraint. Then, in a second step, solve the system of 2 (non-linear...) equations in 2 unknowns $c$ and $d$.

syms  xL xR c d;%symbolic variables declaration
RL=4;RR=1;%radii
a=1;b=1;m=max([RL,RR]);p=2*m+a;
axis([-p+3,p+3,-p-3,p+3]);axis equal;
f=@(x)(b+c*x*x+d*x*x*x);%cubic function
fp=@(x)(2*c*x+3*d*x*x);%its derivative
%The system of 4 equations in 4 unknowns ; sol. in  [XL,XR,C,D]
[XL,XR,C,D]=solve(...
(xL+RL+a)^2+f(xL)^2==RL^2,
(xR-RR-a)^2+f(xR)^2==RR^2,...
xL+f(xL)*fp(xL)==-RL-a,...
xR+f(xR)*fp(xR)==RR+a,...
xL,xR,c,d);% the 4 variables
plot([-p,p],[0,0]);$x axis
plot([0,0],[-p,p]);%y axis
UC=exp(i*(0:0.01:2*pi));%unit circle
plot(-(RL+a)+RL*UC,'k');%circle C_L
plot((RR+a)+RR*UC,'k');%circle C_R
x=-p:0.01:p;%range of values for x
for k=1:length(XL)
   c=C(k);d=D(k);xL=XL(k);xR=XR(k);
   if isreal(c) && isreal(d) && isreal(xL) && isreal(xR)   
   plot(x,b+c*x.*x+d*x.*x.*x,'color',rand(1,3));hold on;
   end;
end;