Newton Raphson method for system of nonlinear equations

Hi, can someone explain what is wrong with the code as the answers are different from the expected answer which I have attached:

x0 = [3; -1.5];  
maxIter = 50; 
tolX = 1e-14; 

x = x0;
xold = x0;
format long
xv = nan(maxIter,1);
yv = xv;
for i = 1:maxIter
    [f,j] = newtraph(x);
    x = x- j\f;
    
    err(:,i) = abs(x-xold);
    xold = x;
    if (err(:,i)<tolX)
        break;
    end
    xv(i,:) = x(1);
    yv(i,:) = x(2);
end

OutputValues = table(xv(~isnan(xv)), yv(~isnan(yv)), 'VariableNames',{'X','Y'})

display(x)
fprintf('Number of iterations: %d \n',i)
fprintf('Final tolerance before accepted tolerance: %d \n' ,err(i))

function  [fval,jac] = newtraph(X)
x = X(1);
y = X(2);

fval(1,1)=x+exp(-x)+y^3;
fval(2,1)=x^2+2*x*y-y^2+tand(x);

jac = [1-exp(-x), 3*y^2;
        2*x+2*y+1+(tand(x))^2, 2*x-2*y];
end

The output of the code is:

OutputValues =

  23×2 table

           X                    Y        
    ________________    _________________

    3.64477577226406    -1.54258686595204
    3.71768338513548    -1.55256145261782
     3.7349599350235    -1.55483810676808
    3.73922377944533     -1.5554091461244
    3.74028579368963    -1.55555181328355
    3.74055090524151    -1.55558745542312
    3.74061712194685    -1.55559635948701
    3.74063366312717    -1.55559858386503
    3.74063779531838    -1.55559913954898
    3.74063882759969    -1.55559927836732
     3.7406390854791    -1.55559931304626
    3.74063914990129    -1.55559932170959
    3.74063916599493    -1.55559932387382
    3.74063917001537    -1.55559932441447
    3.74063917101973    -1.55559932454954
    3.74063917127064    -1.55559932458328
    3.74063917133332    -1.55559932459171
    3.74063917134898    -1.55559932459382
    3.74063917135289    -1.55559932459434
    3.74063917135387    -1.55559932459447
    3.74063917135411    -1.55559932459451
    3.74063917135417    -1.55559932459451
    3.74063917135419    -1.55559932459452


x =

   3.740639171354190
  -1.555599324594516

Number of iterations: 24 
Final tolerance before accepted tolerance: 8.663327e-09 

The steps is the iteration number, the delta(xk-1) is the tolerance which should be less than

tolX = 1e-14;

and the xk and yk are of course the x and y values.

The code works for a previous example that I tried, where the function had different fval and jac equations, but not for this instance.

Answer 1

[Ayush]

Hi,

It seems that your output tolerance is not under the specified "tolX". To improve this, you can use a "norm" (like the Euclidean norm) function to get a single value representing the "distance" between your old and new "x" values and then compare this against your tolerance. Currently, the tolerance check is performed against the vector "err(:,i)" having absolute error calculation directly; instead of it, update the error using the "norm" of "deltaX" Refer to the modified code:

x0 = [3; -1.5];  
maxIter = 50; 
tolX = 1e-14; 

x = x0;
xold = x0;
format long
xv = nan(maxIter,1);
yv = xv;
err = zeros(1, maxIter); % Initialize the error array

for i = 1:maxIter
    [f,j] = newtraph(x);
    deltaX = j\f; % Calculate the change in x
    x = x - deltaX; % Update x
    
    err(i) = norm(deltaX); % Update the error using the norm of deltaX
    if err(i) < tolX
        break;
    end
    xv(i) = x(1); % Store the x and y values
    yv(i) = x(2);
end

% Remove NaN values from xv and yv
xv = xv(~isnan(xv));
yv = yv(~isnan(yv));

OutputValues = table(xv, yv, 'VariableNames',{'X','Y'})

OutputValues = 23x2 table
           X                    Y        
    ________________    _________________

    3.64477577226406    -1.54258686595204
    3.71768338513548    -1.55256145261782
     3.7349599350235    -1.55483810676808
    3.73922377944533     -1.5554091461244
    3.74028579368963    -1.55555181328355
    3.74055090524151    -1.55558745542312
    3.74061712194685    -1.55559635948701
    3.74063366312717    -1.55559858386503
    3.74063779531838    -1.55559913954898
    3.74063882759969    -1.55559927836732
     3.7406390854791    -1.55559931304626
    3.74063914990129    -1.55559932170959
    3.74063916599493    -1.55559932387382
    3.74063917001537    -1.55559932441447
    3.74063917101973    -1.55559932454954
    3.74063917127064    -1.55559932458328

display(x)

x = 2x1

3.740639171354190 -1.555599324594516

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

fprintf('Number of iterations: %d \n',i)

Number of iterations: 24 

fprintf('Final tolerance before accepted tolerance: %d \n' ,err(i)) % Adjusted the print format to show more precision

Final tolerance before accepted tolerance: 3.722328e-15 

function  [fval,jac] = newtraph(X)
x = X(1);
y = X(2);
fval(1,1)=x+exp(-x)+y^3;
fval(2,1)=x^2+2*x*y-y^2+tand(x);
jac = [1-exp(-x), 3*y^2;
        2*x+2*y+1+(tand(x))^2, 2*x-2*y];
end

This resulted in the error of "3.722328e-15", which is less than your defined tolerance. For more information on the "norm" function, refer to the below documentation:

https://www.mathworks.com/help/matlab/ref/norm.html