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.
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".
Using the "trick" is
$$G(\varphi) = F(\varphi) \sin^{2}\varphi.$$
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).