# Transcendental equation: finding root

The following transcendental equation:

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

with constants $X$, $Y$, $R$, $\lambda$, $\varphi_{0}$ represents a solution of the enginnering problem, where:

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

Its derivative is

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

For:

$X=-6.3490e3$, $Y=-1.3663e4$, $R=6378$, $\lambda=-3/4\pi$, $\varphi_{0} = \pi/12$, the $F(\varphi)$ course at interval $\varphi\in[-\pi/2,\pi/2]$ is

Taking into account the above-mentioned facts, which numeric method is appropriate to find the root $\varphi_{0}=-\pi/4$, $F(\varphi_{0})=7.4505e-8$? Unfortunately, the Newton method cannot be applied due to the non-convexity of $F(\varphi)$ (I can be wrong).

Currently, I am using the differential evolution with the objective function

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

which is unnecessary and computationally expensive.

Matlab code:

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)


....................................UPDATED QUESTION....................................

Using the Newton's method:

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


the result depends on the initialization. Unfortunately, for lat = 0.1, the Newton's method converges to $\varphi=9.248346$ (see the last column)

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


However, for lat = -0.1, the root $\varphi=-0.785398$ is found:

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


Hence, the Newton's method is sensitive to the initial guess (positive or negative) that may not be known a priori.

@ Claude Leibovici:

The modified version of $F$ has the form of

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

its course is "flipped".

Using the "trick" is

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

• Why can't you use Newton's method? What other methods can't you use? – Robert Israel Aug 11 '17 at 22:57
• @ Robert: In my opinion, F is not convex at [-pi/2, pi/2] and F'(0) does not exist. – justik Aug 11 '17 at 23:04
• Convexity is not a requirement for Newton's method. The fact that $F$ is undefined at $0$ is not a problem either: just make sure your iterations stay in a sufficiently small interval that contains a root. – Robert Israel Aug 11 '17 at 23:40
• Newton's method with initial point $\varphi = -1$ should work quite well. – Robert Israel Aug 11 '17 at 23:47
• @ Robert: Yes, but only for $\varphi<0$. However, for $\varphi>0$, the initial guess $\varphi = -1$ will not work. See the sample code of the updated question,please. It is very sensitive to the sign of the initial guess which is not known a priori. – justik Aug 12 '17 at 0:08

The problem comes from the discontinuities in particular at $\phi=0$. In fact, your equation has an infinite number of roots and, depending where you start iterating, you could cross many of these discontinuities.

Starting using $\phi_0=0.1$, you converge to the third positive root of the equation.

In order to overcome this problem, I should instead solve $$G(\phi)=F(\phi) \, \sin^2(\phi)=0$$ ignoring, for sure, the trivial solution $\phi=0$. This removes all the discontinuities.

Using your numbers, the plot is quite nice. However, in the range $-\frac \pi 2 \lt \phi \lt \frac \pi 2$ the function $G(\phi)$ goes through a minimum and you need to start on the left of it.

What is amazing is that, if you make a Taylor expansion of $G(\phi)$ up to second order around $\phi=-\frac \pi 4$, solving the quadratic, you find $\phi\approx -0.785398$ which is ... your solution !

Making a Taylor expansion of $G(\phi)$ up to second order around $\phi=0$, solving the quadratic, you find $\phi\approx -0.576201$ which is a very good starting point.

Edit

In any manner, when you know that the solution of $f(x)=0$ is unique and such that $a < x < b$, I strongly suggest the use of combination of bisection and Newton steps.

In the book "Numerical Recipes" (google for it), you could find subroutine rtsafe (Fortran source code here).

I have been using it millions of times (and even improved it implementing Halley and Householder methods).

• @ Claude: Thank you very much for your advice and explanation.This is a simplified example, I try to adapt the proposed solution for the "full" version of the problem. – justik Aug 12 '17 at 11:46
• @justik. You are very welcome ! It is very easy to mimic rtsafe in any computing environment. The trick of removing discontinuities is one of my favoured. I spent most of my life with that ! Please, let me know how it works for the full version. Cheers. – Claude Leibovici Aug 12 '17 at 11:58
• @ Claude: NR is an excellent book, I found the C++ version of the above-mentioned code, which is easier to understand. Removing discontinuities represents a nice trick, it was successful as well as for the entire problem. Merci et Bon soir. – justik Aug 12 '17 at 22:35
• @justik. May I suggest you add a plot of the "new" function in your post ? – Claude Leibovici Aug 13 '17 at 2:24
• @ Claude: The function has been attached. Thanks again. – justik Aug 13 '17 at 22:50