4
$\begingroup$

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

enter image description here

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.

enter image description here

Thanks for your help.

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".

enter image description here

Using the "trick" is

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

enter image description here

$\endgroup$
  • $\begingroup$ Why can't you use Newton's method? What other methods can't you use? $\endgroup$ – Robert Israel Aug 11 '17 at 22:57
  • $\begingroup$ @ Robert: In my opinion, F is not convex at [-pi/2, pi/2] and F'(0) does not exist. $\endgroup$ – justik Aug 11 '17 at 23:04
  • $\begingroup$ 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. $\endgroup$ – Robert Israel Aug 11 '17 at 23:40
  • $\begingroup$ Newton's method with initial point $\varphi = -1$ should work quite well. $\endgroup$ – Robert Israel Aug 11 '17 at 23:47
  • $\begingroup$ @ 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. $\endgroup$ – justik Aug 12 '17 at 0:08
4
$\begingroup$

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).

$\endgroup$
  • $\begingroup$ @ 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. $\endgroup$ – justik Aug 12 '17 at 11:46
  • $\begingroup$ @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. $\endgroup$ – Claude Leibovici Aug 12 '17 at 11:58
  • $\begingroup$ @ 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. $\endgroup$ – justik Aug 12 '17 at 22:35
  • $\begingroup$ @justik. May I suggest you add a plot of the "new" function in your post ? $\endgroup$ – Claude Leibovici Aug 13 '17 at 2:24
  • $\begingroup$ @ Claude: The function has been attached. Thanks again. $\endgroup$ – justik Aug 13 '17 at 22:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.