Multivariate Newton-Raphson method and FindRoot module [closed]

Let's suppose that we have the following equation

Clear["Global*"];

m = 1/2;
V = m/Sqrt[(x - m)^2 + (y - m)^2 + (z - m)^2] + m/Sqrt[(x + m)^2 + (y + m)^2 + (z + m)^2]
+ 1/2*(x^2 + y^2);

Vx = D[V, x];
Vy = D[V, y];
Vz = D[V, z];


Then we can use the FindRoot module for finding a solution of the system $$Vx = Vy = Vz = 0$$ with initial conditions $$(x_0,y_0,z_0) = (1,1,1)$$.

sol = Module[{s = 0, e = 0}, {FindRoot[{Vx == 0, Vy == 0, Vz == 0},
{{x, 1}, {y, 1}, {z, 1}}, WorkingPrecision -> 16,
StepMonitor :> s++, EvaluationMonitor :> e++],
"Steps" -> s, "Evaluations" -> e}]


The results is

{{x -> 0.9466454951651958, y -> 0.9466454951651958, z -> 0.4788183020073119}, "Steps" -> 7, "Evaluations" -> 9}

Now I want to obtain the same result but using my own version of the numerical Newton-Raphson iterator, according to these notes.

Vxx = D[Vx, x];
Vxy = D[Vx, y];
Vxz = D[Vx, z];

Vyx = D[Vy, x];
Vyy = D[Vy, y];
Vyz = D[Vz, z];

Vzx = D[Vz, x];
Vzy = D[Vz, y];
Vzz = D[Vz, z];

ff = Vx;
gg = Vy;
hh = Vz;

fx = Vxx;
fy = Vxy;
fz = Vxz;
gx = Vyx;
gy = Vyy;
gz = Vyz;
hx = Vzx;
hy = Vzy;
hz = Vzz;

J11 = gy*hz - hy*gz;
J12 = hx*gz - gx*hz;
J13 = gx*hy - hx*gy;
J21 = hy*fz - fy*hz;
J22 = fx*hz - hx*fz;
J23 = hx*fy - fx*hy;
J31 = fx*gz - gy*fz;
J32 = gx*fz - fx*gz;
J33 = fx*gy - gx*fy;

det = fx*J11 + fy*J12 + fz*J13;

Nx = ff*J11 + gg*J21 + hh*J31;
Ny = ff*J12 + gg*J22 + hh*J32;
Nz = ff*J13 + gg*J23 + hh*J33;

x0 = 1.;
y0 = 1.;
z0 = 1.;

tol = 10^-15;
n = 20;
iter0 = 0;

Do[
iter = iter0 + i;
Nx0 = Nx /. {x -> x0, y -> y0, z -> z0};
Ny0 = Ny /. {x -> x0, y -> y0, z -> z0};
Nz0 = Nz /. {x -> x0, y -> y0, z -> z0};
det0 = det /. {x -> x0, y -> y0, z -> z0};
x1 = x0 - Nx0/det0;
y1 = y0 - Ny0/det0;
z1 = z0 - Nz0/det0;
dx = Abs[x1 - x0];
dy = Abs[y1 - y0];
dz = Abs[z1 - z0];
Print["i = ", iter, ",  x1 = ", NumberForm[x1, 16], ",  y1 = ",
NumberForm[y1, 16], ",  z1 = ", NumberForm[z1, 16]];
If[Abs[x1 - x0] < tol && Abs[y1 - y0] < tol &&
Abs[z1 - z0] < tol, {Print["The method converges"], Exit[]}];
x0 = x1;
y0 = y1;
z0 = z1;
, {i, 1, n}
];


As you can see, the iterative method diverges and does not provide the results given by FindRoot.

So my question is: Is there any mistake in the implementation of the iterative scheme? Why it fails to provide the same results as FindRoot?

NOTE: For other systems of equations the result of the iterator coincides with that of FindRoot.

closed as off-topic by JimB, m_goldberg, Michael E2, José Antonio Díaz Navas, bbgodfreyJan 22 at 1:34

This question appears to be off-topic. The users who voted to close gave this specific reason:

• "The question is out of scope for this site. The answer to this question requires either advice from Wolfram support or the services of a professional consultant." – m_goldberg, bbgodfrey
If this question can be reworded to fit the rules in the help center, please edit the question.

• I'm voting to close this question as off-topic because while the question is certainly important to the OP, it is unlikely to be of value to someone else. My understanding is that this site is not meant for debugging code that attempts to duplicate existing Mathematica functions. – JimB Jan 16 at 22:23
• Even if your implementation of the inverse of the Jacobian is correct, the method may be still divergent for not-so-nice initial points. Usually, adding a line search can enlarge of the set of "good" initial points tremendously. Btw.: Why don't you use vectors and matrices (and functions instead of /.)? – Henrik Schumacher Jan 17 at 9:01
• @HenrikSchumacher Next you'll be asking why not use Det[], or what's the point of the 10th through 21st lines of code in the second code-block, or Print[], or Exit[] of all things (inside a list no less).... :) It's clearly likely a line-by-line translation of C-like code. Perhaps the OP wants it corrected before refactoring. – Michael E2 Jan 17 at 13:08
• Try x0 = 1.; y0 = 1.; z0 = 1./2; The original starting point diverges, which, as is well known, can happen with Newton's method. – Michael E2 Jan 18 at 0:25
• Then your question is really about FindRoot, which, according to the tutorial, does not use a pure Newton's method, but a damped one with step control. Try FindRoot with Method -> {"Newton", "StepControl" -> None}. – Michael E2 Jan 18 at 14:25

With

V = m/Sqrt[(x - m)^2 + (y - m)^2 + (z - m)^2] +
m/Sqrt[(x + m)^2 + (y + m)^2 + (z + m)^2] + 1/2*(x^2 + y^2);


The system

$$V_x = 0\\ V_y = 0\\ V_z = 0$$

has notoriously three solutions which are for $$m = \frac 13$$

$$X_a\left(\frac 13\right) = \{0.751246,0.751246,0.304505\}\ \ \ \mbox{or}\ \ \ X_b=\{0,0,0\}$$

For $$X_0 > m\{1,1,1\}$$ the iterative procedure converges to $$X_a(m)$$ If $$X_0 < -m\{1,1,1\}$$ converges to $$-X_a(m)$$ otherwise converges to $$X_b$$

Follows a very basic script which shows this fact.

X = {x, y, z};
alpha = m + 0.1;
minerror = 10^-10;
iH = Inverse[H];
X0 = alpha {1., 1., 1.};

For[i = 1, i <= 20, i++,
iH0 = iH /. Thread[X -> X0];
f0 = f /. Thread[X -> X0];
error = Max[Abs[f0]];
Print[i, " " , X0, " ", error];
If[error < minerror, Break[]];
X1 = X0 - iH0.f0;
X0 = X1
]
`

Resuming, the iterative process converges to $$X_a(m)$$ when $$X_0 > m\{1,1,1\}$$. If $$X_0 < -m\{1,1,1\}$$ converges to $$-X_a(m)$$, otherwise, converges to $$X_b$$

NOTE

The iterative process can be described as

$$X_{k+1} = X_k - H_k^{-1}\cdot \nabla V_k$$

where $$X = \{x,y,z\},\ \ H = \nabla^2 V$$

Attached the three solutions for $$m=\frac 12$$. In light red $$V_x = 0$$ in light green $$V_y = 0$$ and in light yellow $$V_z = 0$$

• This does not solve my problem. I want the solution of the system for m =1/2 not 1/3. Moreover, I am asking why the custom iterative scheme diverges.. – Vaggelis_Z Jan 16 at 21:59
• @Vaggelis_Z Please. Substitute in the script m=1/2 and it will converge to $X_a(1/2) = \{0.946645,0.946645,0.478818\}$ – Cesareo Jan 16 at 22:02
• OK but again I need to know why my custom iterative scheme does not converge to the same root. Any ideas? – Vaggelis_Z Jan 16 at 22:07
• @Vaggelis_Z Having MATHEMATICA with all it's symbolic processing power, the most indicated practice is to use this power to avoid transcription errors that may arise when manipulating formulas obtained with intermediate steps as could be the case with your script. – Cesareo Jan 16 at 22:15

The iterates of Newton's method for the equation $$F(x) = 0$$ with a initial guess $$x_0$$ look like $$x_{k+1} = x_k - DF(x_k)^{-1} \, F(x_k).$$ If $$x_0$$ is chosen badly, then there may happen large jumps $$x_{k+1} - x_k$$, so large that you can easily miss the basin of attraction.

One can remedy this by a picking an appropriate step size $$0 < \tau_k \leq 1$$ and use the iterates $$x_{k+1} = x_k - \tau_k \, DF(x_k)^{-1} \, F(x_k).$$

Choosing a fixed $$\tau_k = q$$ with $$0< q< 1$$ is often referred to as damped Newton method. Alas, this slows down the convergence rate from $$q$$-quadratic convergence to only $$q$$-linear convergence (if I recall correctly).

Line search methods try to automize the choice of $$\tau_k$$ in such a way that one obtains $$\tau_k \approx 1$$ if $$x_k$$ is close to a solution in order to obtain optimal convergence rate.

Armijo line search (also called backtracking line search) applies the following strategy: It beliefs that the search direction $$u_k = - DF(x_k)^{-1} \, F(x_k)$$ is already good and picks $$\tau_k$$ such that a certain merit function $$\varphi \colon [0,\infty[ \to \mathbb{R}$$ is sufficiently reduced in the sense that $$\varphi(\tau_k) < \varphi(0) - \sigma \varphi'(0) \, \tau_k$$ with some parameter $$\sigma \in ]0,1[$$.

If one applies Newton's method for minimization, i.e., if one tries to solve the equation $$Df(x)=0$$ with a given objective function $$f$$, then one usually choses the merit function as $$\varphi(t) = f(x_k + t \, u_k)$$. Then one starts with an initial guess $$t$$ for $$\tau_k$$, say $$t = 1$$, checks whether $$\varphi(t) < \varphi(0) - \sigma \varphi'(0) \, t$$ is fulfilled: If it is fulfilled, one choose $$\tau_k = 1$$; otherwise, one scales $$t$$ down (i.e. $$t \leftarrow \gamma \, t$$ with some $$0<\gamma<1$$) and tries again... until one finds a $$t$$ satisfying the inequality. If $$f$$ is sufficiently smooth, existence of such a $$t$$ is guaranteed.

For general equations $$F(x) = 0$$, this won't work since there is no function $$f$$. Then one may use the merit function $$\varphi(t) = \|F(x_k + t\, u_k)\|^2$$. Notice that $$\varphi'(0) = 2 \langle F(x_k) , DF(x_k) \, u_k \rangle = -2 \|F(x_k)\|^2$$, so that we obtain the following as condition of sufficient reduce of the merit function: $$\|F(x_k + t \, u_k)\|^2 < \|F(x_k)\|^2 (1- 2 \, \sigma \, t).$$ In fact, this can be satisfied only if $$\sigma < \frac{1}{2}$$. In practice, $$\sigma$$ can be chosen rather small; often $$\sigma = \frac{1}{1000}$$ and $$\gamma = \frac{1}{4}$$ work quite fine. This basically ensures $$\|F(x_{k+1})\|^2 < \|F(x_k)\|^2$$ which is fairly reasonable because we want to minimize $$\|F(x)\|^2$$. Since we want maximal convergence rate when $$x_k$$ is close to a solution, we should always use $$t = 1$$ as initial guess for $$\tau_k$$.