Ecuación trascendental: encontrar la raíz

Question

La siguiente ecuación trascendental:

$$F(\varphi)=(X^{2}+(S+\rho-Y)^{2}-\varrho^{2} = 0,$$

con las constantes $X$ , $Y$ , $R$ , $\lambda$ , $\varphi_{0}$ representa una solución del problema de enginnering, donde:

$$\rho =\frac{R}{\tan\varphi}, $$ $$\varepsilon=\lambda\sin\varphi,$$ $$S = R (\varphi - \varphi_{0}).$$

Su derivado es

$$F^{\prime}(\varphi)=2\varrho(R-\frac{S}{\tan\varphi}+\frac{Y}{\tan\varphi}).$$

Para:

$X=-6.3490e3$ , $Y=-1.3663e4$ , $R=6378$ , $\lambda=-3/4\pi$ , $\varphi_{0} = \pi/12$ El $F(\varphi)$ curso en el intervalo $\varphi\in[-\pi/2,\pi/2]$ es

Teniendo en cuenta los hechos mencionados, ¿qué método numérico es el adecuado para encontrar la raíz $\varphi_{0}=-\pi/4$ , $F(\varphi_{0})=7.4505e-8$ ? Lamentablemente, el método de Newton no puede aplicarse debido a la no convexidad de $F(\varphi)$ (Puedo equivocarme).

Actualmente, estoy utilizando la evolución diferencial con la función objetivo

$$\phi = \left|F(\varphi)\right|,$$

lo cual es innecesario y costoso desde el punto de vista computacional.

Gracias por su ayuda.

Código Matlab:

X = -6.349068367938568e+03
Y = -1.366383308370021e+04
R = 6378 
lam = -3/4*pi
phi0 = pi/12;

phi = -pi/2:0.01:pi/2
S = R * (phi - phi0);
rho = R ./ tan (phi);
F = (X.^2 + (S + rho - Y).^2 - rho.^2);
plot (phi, F)

....................................PREGUNTA ACTUALIZADA....................................

Utilizando el método de Newton:

phi = -0.1
for i = 1:10
    S = R * (phi - phi0);
    rho = R / tan (phi);
    F = (X^2 + (S + rho - Y)^2 - rho^2);
    dF = 2 * rho * (R - S/tan(phi) + Y / tan(phi));
    phi = phi - F / dF;

    fprintf('%d  %f  %f  %f  %f %f \n',i, S, rho, F, dF, phi );

end

el resultado depende de la inicialización. Por desgracia, para lat = 0,1, el método de Newton converge a $\varphi=9.248346$ (véase la última columna)

1  -1031.956495  63567.258132  1805822494.834060  -15195043945.052212   0.218843 
2  -273.976728  28677.447286  987572721.972760  -3087236602.757069   0.538732 
3  1766.274678  10670.780257  607701473.229591  -414826683.790066   2.003684 
4  11109.743042  -2947.409542  508004993.852166  -105083383.937407   6.837988 
5  41942.932166  10291.527860  4276980158.270769  -1715576361.875511   9.331016 
6  57843.465687  -67824.110109  -4546233364.176280  -104013936332.721650   9.287308 
7  57564.696539  -46103.099453  -1453897868.747859  -48062515428.430321   9.257058 
8  57371.761140  -37670.492327  -265525298.526285  -32090495323.084854   9.248784 
9  57318.987870  -35864.942274  -12718012.164382  -29088657418.251171   9.248347 
10  57316.199310  -35774.200614  -32169.163826  -28941685056.659229   9.248346

Sin embargo, para lat = -0,1, la raíz $\varphi=-0.785398$ se encuentra:

1  -2307.556495  -63567.258132  -1274499043.514500  -15200453311.959137   -0.183846 
2  -2842.327053  -34300.316046  -584946491.969031  -4429892388.511618   -0.315891 
3  -3684.512009  -19514.385761  -249583123.903867  -1440596169.114301   -0.489141 
4  -4789.499790  -11982.287733  -93605169.378424  -552386466.495069   -0.658597 
5  -5870.289792  -8241.804207  -27415729.600543  -271138926.206457   -0.759710 
6  -6515.189903  -6714.387159  -4583747.486028  -186709319.718563   -0.784261 
7  -6671.770953  -6392.526587  -194284.152926  -171140559.538748   -0.785396 
8  -6679.011460  -6378.029043  -388.102944  -170457334.568212   -0.785398 
9  -6679.025981  -6378.000000  -0.001556  -170455967.400742   -0.785398 
10  -6679.025982  -6378.000000  0.000000  -170455967.395259   -0.785398 

Por lo tanto, el método de Newton es sensible a la conjetura inicial (positiva o negativa) que puede no ser conocida a priori.

---

@ Claude Leibovici:

La versión modificada de $F$ tiene la forma de

$$F(\varphi)=[(X^{2}+(S+\rho-Y)^{2}-\varrho^{2} + 1/ \varphi^{2}]\cos2\varphi = 0.$$

su curso está "volteado".

El uso del "truco" es

$$G(\varphi) = F(\varphi) \sin^{2}\varphi.$$

Answer 1

El problema viene de las discontinuidades en particular en $\phi=0$ . De hecho, su ecuación tiene un número infinito de raíces y, dependiendo de dónde empiece a iterar, podría cruzar muchas de estas discontinuidades.

Empezar a usar $\phi_0=0.1$ , se converge a la tercera raíz positiva de la ecuación.

Para superar este problema, debería resolver en cambio $$G(\phi)=F(\phi) \, \sin^2(\phi)=0$$ ignorando, por supuesto, la solución trivial $\phi=0$ . Esto elimina todas las discontinuidades.

Usando tus números, la trama es bastante bonita. Sin embargo, en el rango $-\frac \pi 2 \lt \phi \lt \frac \pi 2$ la función $G(\phi)$ pasa por un mínimo y hay que empezar por la izquierda del mismo.

Lo que es sorprendente es que, si se hace una ampliación de Taylor de $G(\phi)$ hasta el segundo orden en torno a $\phi=-\frac \pi 4$ resolviendo la cuadrática, se encuentra $\phi\approx -0.785398$ que es... ¡su solución!

Hacer una ampliación de Taylor de $G(\phi)$ hasta el segundo orden en torno a $\phi=0$ resolviendo la cuadrática, se encuentra $\phi\approx -0.576201$ que es un muy buen punto de partida.

Editar

De cualquier manera, cuando se sabe que la solución de $f(x)=0$ es único y tal que $a < x < b$ Sugiero encarecidamente el uso de la combinación de la bisección y los pasos de Newton.

En el libro "Numerical Recipes" (búsquelo en Google), puede encontrar la subrutina rtsafe (Código fuente Fortran aquí ).

Lo he utilizado millones de veces (e incluso lo he mejorado aplicando los métodos Halley y Householder).