## Welcome to the World of Modelling and Simulation

### What is Modelling?

This blog is all about system dynamics modelling and simulation applied in the engineering field, especially mechanical, electrical, and ... ### 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.