problem with programm that use Newton-Raphson method to solve nonlinear equation system

I have to solve a nonlinear equation system with the Newton-Raphson method, and I wrote a code according to a literature. I have tried this code with simply equations, and it works very well. But wenn I use the same code to solve my own equations, here is the problem: the results changed with different precisions(for example 1e-6 and 1e-10) or guess values, and the diefference between the results is very big (1e-3 times). How can I fix that?

The function to creat the equation system (fun=myfun(g)) is very complicated and very long, I will paste it if it's necessary.

% Jacobi matrix
function dfun=Jacobi(g);
syms U_mn V_mn W_mn phi_xmn phi_ymn
g=[U_mn; V_mn; W_mn; phi_xmn; phi_ymn];
fun=myfun(g);
dfun=[diff(fun,U_mn);diff(fun,V_mn);diff(fun,W_mn);diff(fun,phi_xmn);diff(fun,phi_ymn)];
dfun=conj(dfun');

% Newton-Raphson method
clc
clear all
g0=[0.000001 0.000001 0.000001 0.000001 0.000001]; % guess value
eps=1e-6; % precision
it_max=50; % maximal iteration steps
con=0; % convergence
iter=1
for i=1:it_max;
  fun=myfun(g0);
  dfun=Jacobi(g0);
  fun=subs(fun,{'U_mn' 'V_mn' 'W_mn' 'phi_xmn' 'phi_ymn'},{g0(1) g0(2) g0(3) g0(4) g0(5)});
  dfun=subs(dfun,{'U_mn' 'V_mn' 'W_mn' 'phi_xmn' 'phi_ymn'},{g0(1) g0(2) g0(3) g0(4) g0(5)});
  g=g0-fun/dfun;
  if norm(g-g0)<eps
      con=1;
      break;
  end
  g0=double(g)
  iter=iter+1
end
format long
g=double(g)

Answer 1

[Riccardo Scorretti]

Hi Mi,

it is difficoult to give a precise answer without a deeper knowledge of the equation which has to be solved. Unfortunately, convergency issues are very often found with Newton-Raphson, or any other method to solve nonlinear equations.

Up to my (very limited) experience, I can give you a few advices:

1. You could try to use Matlab fsolve function instead of writing a custom nonlinear solver.
2. From your code, I see that the only convergency test is on the absolute variation of the candidate solution (abs(g-g0)<eps). Perhaps you ought to test also the residual. Kelley suggests to test , where x = current iterate, = initial quest, the equation to be solved is , and and are respectivel the relative and absolute tolerances.
3. Most importantly, you could implement strategies like Armijo rule, also known as "backtracking". In practice, at each iteration compute the tentative direction (in your code -fin/dfun). Then, instead of just computing the next iterate g=g0-fun/dfun, an inner loop is added so as to ensure that the residual is reduced at the next iteration. This can greatly improve the convergency properties of Newton-Raphson method.

I can provide the following code, which is extremely simplistic and not at all optimized, which illustrate these concepts on a simple case. The (scalar) equation which is solved is: , the solution of which is obviously . It can be observed that the method diverges is Armijo rule is not used.

% Equation to solve
% -----------------
f = @(x) atan(x) - 1.0;
df = @(x) 1.0/(x^2+1.0);

% Plot the graphic
figure
x = linspace(-5, 5, 1024);
plot(x, f(x)) ; grid on ; hold on
xlabel('x') ; ylabel('f(x)');

% Set the parameters of the solver
tau_a = 1.0e-8;     % = absolute tolerance
tau_r = 1.0e-8;     % = relative tolerance
x = 40.0;           % = initial guess
maxit = 40;         % = max. number of iteraionts
use_armijo = true;  % = use Armijio rule

f0 = f(x);
for n = 1 : maxit
  delta_x = - f(x)/df(x);
  
  % Armijio rule (it improves the global convergency)
  while use_armijo  &&  abs(f(x)) < abs(f(x+delta_x))
    delta_x = 0.7 * delta_x;
  end
  
  x = x + delta_x;
  fprintf('x = %20.15f\n', x);
  
  % Exit test
  if abs(f(x)) < tau_r*(abs(f0)) + tau_a
    disp('*** Convergence ***');
    plot(x, f(x), 'r+', 'MarkerSize', 12);
    break;
  end
end

x =    4.737878027780923
x =    1.820188841802746
x =    1.525092387763806
x =    1.556934889460351
x =    1.557407623013678
x =    1.557407724654897

*** Convergence ***

I may suggest to read the book: C.T. Kelley, Solving nonlinear equations with Newton's method, SIAM (ISBN 0-89871-546-6). It contains many Matlab codes, and it is only 100 pages long.

I hope it helps.