Polinomio cúbico que une suavemente dos círculos

Question

Dados dos círculos con radios $R_L$ y $R_R$ y centros en $(-(R_L+a),\,0)$ y $(R_R+a,\,0)$ respectivamente, encontrar un polinomio cúbico $p(x)=b+cx^2+dx^3$ que conecta suavemente los dos círculos.

$b$ es un parámetro por lo que $p(0)=b$ y el término lineal del polinomio se omite porque queremos $\frac{\mathrm{d}p}{\mathrm{d}x}\Big|_{x=0}=0$ .

Mi intento de solución. Deja que $\mathrm{C}_{L,R}$ las ecuaciones de la mitad superior del $L,R$ círculos. Formulo dos ecuaciones que relacionan $\mathrm{C}_{L,R}$ y $p$ y dos ecuaciones que relacionan $\frac{\mathrm{d}}{\mathrm{d}x}\mathrm{C}_{L,R}$ y $\frac{\mathrm{d}p}{\mathrm{d}x}$ . Sea $x_{L,R}$ sean los puntos donde $p(x)$ y $\mathrm{C}_{L,R}(x)$ se cruzan, entonces:

$$ \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 $$ así que tengo un sistema de cuatro ecuaciones no lineales con cuatro incógnitas $(x_L,\,x_R,\,c,\,d)$ . He codificado un simple método de Newton para el sistema y funciona bien para algunas combinaciones de parámetros $(a,\,b,\,R_L,\,R_R)$ cuando la suposición inicial es lo suficientemente cercana, especialmente cuando $|R_L-R_R|$ no es demasiado grande y utilizo un amortiguamiento constante para las iteraciones de Newton. Puedo encontrar conjeturas iniciales que me parecen buenas a través de una interfaz gráfica que he codificado. Sin embargo, como $|R_L-R_R|$ se hace más grande, el solucionador falla espectacularmente para converger incluso con conjeturas iniciales muy cercanas y amortiguación muy pequeña. (Debo añadir que en realidad estoy tomando el cuadrado de las ecuaciones para evitar las raíces cuadradas de los números negativos durante las iteraciones de Newton).

Mi pregunta es triple:

a) ¿qué otro método o modificación puedo utilizar para que el solucionador sea más estable?

b) este problema me parece que debería estar resuelto en algún sitio, ¿conoces alguna referencia?

c) de manera más general, ¿hay alguna razón para que esto falle tan horriblemente como lo hace cuando $|R_L-R_R|>>1$ ?

Answer 1

Aquí hay un método que da todas las soluciones de forma determinista.

Observe la siguiente figura:

Fig. 1 : El conjunto exhaustivo de 8 curvas de solución en el caso $a=1, \ b=1, \ R_L=3, \ R_R=1$ . Tenga en cuenta que algunas tangencias son internas a los círculos.

¿Cómo se ha obtenido este resultado? Utilizando un sistema de álgebra computacional con un sistema bastante simple de 4 ecuaciones polinómicas en las 4 incógnitas $c,d,x_L,x_R$ donde (sus anotaciones) los dos últimos parámetros son las abscisas de los puntos de tangencia de la curva cúbica con ecuación

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

con el círculo izquierdo y derecho respectivamente.

Aquí están:

$$\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}$$

Las dos primeras ecuaciones expresan que $M_L=(x_L,f(x_L))$ y $M_R=(x_R,f(x_R))$ pertenecen a sus respectivos círculos.

Expliquemos ahora la tercera ecuación. Sea $C_L(-(R_L+a),0)$ denota el centro del primer círculo. La pendiente de la tangente en $M_L$ es $f'(x_L)$ . $\vec{C_LM_L}=\binom{x_L+R_L+a}{f(x_L)-0}$ siendo ortogonal a esta tangente, tiene una pendiente $-\dfrac{1}{f'(x_L)}$ (ver allí ). Un razonamiento similar para la cuarta ecuación.

Observaciones :

1) No he encontrado las limitaciones que mencionas en el caso de una gran brecha entre $R_L$ y $R_R$ . Por ejemplo, cuando $R_L=100$ y $R_R=1$ (con $a=1$ y $b=0$ ), se encuentran 14 soluciones...

2) La tercera ecuación se puede escribir : $f(x_L)f'(x_L)+x_L+R_L+a=0$ : esto evita una posible división por cero. Lo mismo para la cuarta ecuación.

Otro caso (Fig. 2), esta vez con $R_L=R_R$ mostrando algunas soluciones espurias bajo la forma de ... curvas de segundo grado, es decir, parábolas (efectivamente se puede tener $d=0$ ) :

Fig. 2.

Otro caso, con cuatro soluciones (una de ellas parece una parábola pero no lo es):

Fig. 3.

Aquí está el programa Matlab que ha dado la figura 1
(tiempo de ejecución : unos 30 segundos en mi ordenador, bastante lento). Tenga en cuenta las condiciones "isrealistas" (queremos que las diferentes incógnitas sean reales) : de hecho en el primer caso, hay $40$ soluciones, $32$ de que sean espurios, es decir, con coeficientes complejos...

Observación final : De hecho, existe una forma diferente de resolver el sistema (1) haciéndolo en dos pasos ; primero, expresando el hecho de que las ecuaciones (i) y (iii) tienen una raíz común $x_L$ dando (usando un "resultante" ) una primera restricción (no lineal) entre $c$ y $d$ . Haciendo lo mismo para las ecuaciones (ii) y (iv) para la raíz común $x_R$ , da lugar a una segunda restricción. A continuación, en un segundo paso, resolver el sistema de 2 ecuaciones (no lineales...) en 2 incógnitas $c$ y $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;