A MATLAB Program to Solve a System of Non Linear Equations by Newton-Raphson Method

Problem Statement: Develop a MATLAB program to solve a system of nonlinear equations by the Newton-Raphson method, and then, test the code with the following equations:

exp(2 x₁) –x₂ -4 =0.0 
x₂ – x₃² -1 =0.0 
x₃ – sin x₁ =0.0 

Solution:
The following program is developed in MATLAB to implement the Newton-Raphson algorithm. A function named "newtonraphson" is created, which applies the technique.

% Newton-Raphson method applied to a system of non-linear equations f(x)=0,
% given the jacobian function J, where J = del(f1,f2,...,fn)/del(x1,x2,...,xn)
% and x=[x1;x2;...;xn], f=[f1;f2;...;fn]; x0 is the initial guess to run
% the program for solution

function [x,iteration] = newtonraphson(x0,f,J)

N = 100; % Maximum number of iterations definition

epsilon = 1e-15; % Tolerance Definition

maxvalue = 10000.0; % Definition of value for divergence

xxx = x0; % Initial guess for the solution

while (N>0)
   
    JJJ = feval(J,xxx);    % "feval" is used here to evaluate functions

    if abs(det(JJJ))<epsilon
        
        error('newtonraphson - Jacobian is singular - consider new x0');
        
        abort;
        
    end;
    
xn = xxx - inv(JJJ)*feval(f,xxx);  % Newton-Raphson algorithm in matrix form

    if abs(feval(f,xn))<epsilon

        x=xn;
        
        iteration = 100-N;

        return;

    end;
    
    if abs(feval(f,xxx))>maxvalue
        
        iteration = 100-N;

        disp(['iterations = ',num2str(iteration)]);
                
        error('Solution diverges');

        abort;

    end;

    N = N - 1;

    xxx = xn;

end;

error('No convergence after 100 iterations.');

abort;

end  % Function end

The following function determines the derivatives:

function [f] = funcdef(x)

% funcdef(x) = 0, with x = [x(1);x(2);x(3)] represents a system of 3 nonlinear equations

f1 = exp(2*x(1))-x(2)-4;

f2 = x(2)-(x(3)*x(3))-1;

f3 = x(3)-sin(x(1));

f = [f1;f2;f3]; % End of function



The following function determines the Jacobian:


function [J] = jacobmatrix(x)

% Evaluates the Jacobian of a 3x3 for a system of 3 nonlinear equations

J(1,1) = 2*exp(2*x(1));
J(1,2) = -1; 
J(1,3) = 0;

J(2,1) = 0; 
J(2,2) = 1; 
J(2,3) = -2*x(3);

J(3,1) = -cos(x(1)); 
J(3,2) = 0; 
J(3,3) = 1;

end  % Function end

Program Outputs:

>> [x,iteration] = newtonraphson([1;1;1],'funcdef','jacobmatrix')
x =
    0.8590
    1.5733
    0.7572
iteration =
     4

>> funcdef([0.85;1.57;0.75])
ans =
   -0.0961
    0.0075
   -0.0013

So, the solutions are 0.85, 1.57 and 0.75.