Problem with convergence using newtons method

Dear All,

I have problem with convergence using Newtons'method. I am trying to apply Newton's method in Matlab for my problem, and I wrote a script like below for my problem:

% Solve the system of non-linear equations.
% x^2 + y^2 = 2z
% x^2 + z^2 =1/3
% x^2 + y^2 + z^2 = 1
% using Newton’s method having tolerance = 10^(−5)

clc;clear;close all;
syms x y z;

fn = [x^2+y^2-2*z ; x^2+z^2-(1/3);x^2+y^2+z^2-1];
FN = matlabFunction(fn);
jacob_fn = [diff(fn(1) , x) diff(fn(1) , y) diff(fn(1) , z);diff(fn(2) , x) diff(fn(2) , y) diff(fn(2) , z);diff(fn(3) , x) diff(fn(3) , y) diff(fn(3) , z)];
jacob_FN = matlabFunction(jacob_fn);
error = 10^-10 ;
xyz0 = [1 ;1 ;0.1] ;
fnxyz0 = feval(FN,xyz0(1),xyz0(2),xyz0(3));
i = 0;
fprintf('      Iteration|    x     |     y      |    z      | Error      | \n')
while true
    norm1 = norm(fnxyz0);
    fprintf('%10d     |%10.4f| %10.4f | %10.4f| %10.4d |\n',i,xyz0(1),xyz0(2),xyz0(3),norm1)
    jacob_fnxyz0 = feval(jacob_FN,xyz0(1),xyz0(2),xyz0(3));
    H = jacob_fnxyz0\fnxyz0;
    xyz0 = xyz0 - H;
    fnxyz0 = feval(FN,xyz0(1),xyz0(2),xyz0(3));
    i = i + 1 ;
    norm1 = norm(fnxyz0);
    if norm(fnxyz0) < error
        fprintf('%10d     |%10.4f| %10.4f | %10.4f| %10.4d |\n',i,xyz0(1),xyz0(2),xyz0(3),norm1)
        break
    end
    
end

The result this very good:

Iteration| x | y | z | Error |

0 | 1.0000| 1.0000 | 0.1000| 2.1721e+00 |

1 | 0.6258| 0.8333 | 0.4591| 4.3429e-01 |

2 | 0.4432| 0.8167 | 0.4149| 6.0297e-02 |

3 | 0.4041| 0.8165 | 0.4142| 2.6538e-03 |

4 | 0.4022| 0.8165 | 0.4142| 6.2251e-06 |

5 | 0.4022| 0.8165 | 0.4142| 3.4577e-11 |

But when I want to use for my complex function, the results was not converge. The code as following:

clc; clear; close all;
theta = deg2rad([70 80]);      
phi0 = [108.25 65.22];
for i=1:2
    if phi0(i)>90
        phi0(i)=phi0(i)-180;
    end
end
phi=deg2rad(phi0);

% Expected result  n=0.31  k=3.428;
% 
iteration = 5;

syms n k;
u1=sqrt(0.5*((n^2-k^2-sin(theta(1))^2)+sqrt((n^2-k^2-sin(theta(1))^2)^2+4*n^2*k^2)));
v1=sqrt(0.5*(-(n^2-k^2-sin(theta(1))^2)+sqrt((n^2-k^2-sin(theta(1))^2)^2+4*n^2*k^2)));
a1=2*v1*cos(theta(1));
b1=u1^2+v1^2-cos(theta(1))^2;
c1=2*v1*cos(theta(1))*(n^2-k^2-2*u1^2);
d1=u1^2+v1^2-(n^2+k^2)^2*cos(theta(1))^2;
phi1 = ((a1*d1-b1*c1)/(a1*c1+b1*d1))- tan(phi(1));

u2=sqrt(0.5*((n^2-k^2-sin(theta(2))^2)+sqrt((n^2-k^2-sin(theta(2))^2)^2+4*n^2*k^2)));
v2=sqrt(0.5*(-(n^2-k^2-sin(theta(2))^2)+sqrt((n^2-k^2-sin(theta(2))^2)^2+4*n^2*k^2)));
a2=2*v2*cos(theta(2));
b2=u2^2+v2^2-cos(theta(2))^2;
c2=2*v2*cos(theta(2))*(n^2-k^2-2*u2^2);
d2=u2^2+v2^2-(n^2+k^2)^2*cos(theta(2))^2;
phi2 = ((a2*d2-b2*c2)/(a2*c2+b2*d2))-tan(phi(2));


phi = [phi1 ; phi2]
PHI = matlabFunction(phi);
jacob_phi = [diff(phi1,n) diff(phi1,k) ; diff(phi2,n) diff(phi2,k) ];
jacob_PHI = matlabFunction(jacob_phi);
error = 10^-2 ;
n0=0.2; k0=3;
nk = [n0 ;k0] ;
PHI1 = feval(PHI, nk(1), nk(2));
i = 0;
fprintf('      Iteration|    n     |     k      | Error      | \n')
norm1 = norm(PHI1);
fprintf('%10d     |%10.4f| %10.4f | %10.4d |\n',i, n0, k0, norm1)

while true
    
    jacob_PHI1 = feval(jacob_PHI, nk(1), nk(2));
    nk = nk - jacob_PHI1\PHI1;
    PHI1 = feval(PHI, nk(1), nk(2));
    i = i + 1 ;
    norm1 = norm(PHI1);
    fprintf('%10d     |%10.4f| %10.4f | %10.4d |\n',i, nk(1), nk(2), norm1)
    
    if norm(PHI1) < error || i == iteration
        break
    end
end

The results as below:

Iteration| n | k | Error |

0 | 0.2000| 3.0000 | 3.0961e+00 |

1 | 3.8837| 6.8593 | 4.4812e+00 |

2 | 8.6859| 82.1368 | 3.7298e+00 |

3 |6838704.4357| 1468648.3375 | 3.7268e+00 |

4 |-169196838750987247987720192.0000| 76185516560718833553244160.0000 | 3.7268e+00 |

Can anyone help me figure out what I'm doing wrong?

Answer 1

[Torsten]

Use

PHI = matlabFunction(phi,'Vars',{n,k})
jacob_PHI = matlabFunction(jacob_phi,'Vars',{n,k})
n0 = -0.3; k0 = 3.4;

instead of

PHI = matlabFunction(phi);
jacob_PHI = matlabFunction(jacob_phi);
n0=0.2; k0=3;

[chitua Bella]

Dear Torsten,

thanks for your comment. But my problem is for convergence of results. Before I found the problemt that you mentioned.

Now the convergence of data when I choose a very near initial data for {n,k} to the expected value the convergence is good but when I a little bit change the initial data the result not convege, Is there any modification for Newton's method. following I present 3 examples:

we expect to get a possitive results for n and k.

n0=0.2; k0=3.3;

Iteration| n | k | Error |

0 | 0.20000| 3.30000 | 4.47821e-01 |

1 | -0.30099| 3.45305 | 7.09084e-02 |

2 | -0.28789| 3.42954 | 1.44146e-03 |

3 | -0.30195| 3.42871 | 8.91148e-05 |

4 | -0.30162| 3.42871 | 4.60331e-08 |

5 | -0.30162| 3.42871 | 1.11819e-14 |

n0=0.2; k0=3.2;

Iteration| n | k | Error |

0 | 0.20000| 3.20000 | 8.52875e-01 |

1 | -1.52569| 3.49958 | 9.24715e-01 |

2 | -0.14635| 3.61376 | 4.76446e-01 |

3 | 1.14539| 3.47517 | 5.77642e-01 |

4 | 0.33225| 3.48685 | 1.74102e-01 |

5 | 0.22426| 3.43395 | 5.45082e-03 |

6 | 0.31504| 3.42871 | 3.72147e-03 |

7 | 0.30188| 3.42871 | 7.41560e-05 |

8 | 0.30162| 3.42871 | 2.81988e-08 |

9 | 0.30162| 3.42871 | 6.28037e-15 |

n0=0.2; k0=3.1;

Iteration| n | k | Error |

0 | 0.20000| 3.10000 | 1.35968e+00 |

1 | -3.35338| 3.56686 | 4.02760e+00 |

2 | 2.99291| 7.41712 | 5.22637e+00 |

3 |-130.84845| 88.46015 | 3.73519e+00 |

4 |-7886902.27161| 19696183.10014 | 3.72679e+00 |

5 |6907498755769073093591433216.00000| 6588295186410552391666499584.00000 | 3.72679e+00

[Torsten]

Is there any modification for Newton's method ?

Newton's method is the standard method to determine the zeros of systems of nonlinear equations. It has a certain limited region of attraction. If your initial guesses lie outside: bad luck.