## 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 ... ### Explicit Euler Method to Solve System of ODEs in MATLAB

In this tutorial, I am going to show a simple way to solve system of first order ordinary differential equations (ODE) by using explicit Euler method. Let' say, we have following three first order ODEs.

Here, we have 3 ODEs, 3 dependent variables (x, y, z), and 1 independent variable, t. The initial conditions when time is zero are, x(0) = y(0) = z(0) =1. Our goal is to solve these differential equations with explicit Euler approach and plot the solutions afterwards.

Step 1: Define the Equations
The first step is to define all the differential equations in MATLAB. I did this by using MATLAB function handle, which is shown below.

Step 2: Choose a Numerical Approach
The next step is to select a numerical method to solve the differential equations. In this example, we will use explicit Euler method. I have created a function to implement the algorithm. The following image shows the application of the explicit Euler method.

Step 3: Call the Function
This is the final step where we need to call the function. In the previous step, we have created a function, which we are going to call in this step to solve the equations.

So, that's it. We are done with the process, and we can now visualize the solution. I have also attached the MATLAB codes for this problem at the end.

MATLAB Program:

close all;

clc;
format long;

% Defining the three differential equations of the problem
f = @(t,y) [-5*y(1)+5*y(2); 14*y(1)-2*y(2)-y(1)*y(3); -3*y(3)+y(1)*y(2)];
[x,y] = explicit_euler(f,[0,5],[1;1;1],0.001); % Calling Euler function
plot(x,y);
title('When Time Step is 0.001');
legend('x(t)', 'y(t)', 'z(t)', 'Location', 'NorthEast')
xlabel('t')
ylabel('Solutions')

function [x, y] = explicit_euler( f, xRange, y_initial, h )
% This function uses Euler’s explicit method to solve the ODE
% dv/dt=f(t,v); x refers to independent and y refers to dependent variables
% f defines the differential equation of the problem
% xRange = [x1, x2] where the solution is sought on
% y_initial = column vector of initial values for y at x1
% numSteps = number of equally-sized steps to take from x1 to x2
% x = row vector of values of x
% y = matrix whose k-th column is the approximate solution at x(k)
x(1) = xRange(1);
numSteps = ( xRange(2) - xRange(1) ) /h ;
y(:,1) = y_initial(:);
for k = 1 : numSteps
x(k + 1) = x(k) + h;
y(:,k+1) = y(:,k) + h * f( x(k), y(:,k) );
end

#ExplicitEulerMethod #Matlab #NumericalMethod #Blog #Blogger

### Discussion on Numerical Integration Approaches with MATLAB

In this tutorial, we are going to implement four numerical integration schemes in MATLAB. The methods are:

i)   Midpoint Rectangle Rule
ii)  Trapezoidal Rule
iii) Simpson’s 1/3 Rule
iv) Simpson’s 3/8 Rule

Let's say, we would like to integrate the following error function where the upper interval, a is 1.5. At first, we will evaluate the true integral of the error function, and then, we will find the integral numerically with the above methods (i) -(iv).

Finally, we will compare each method with the true value and make some comments based on the deviations if we have any. For each method, we will find out the following parameters:

Here, n refers to the numbers, 1,2,4,8,16, ...,128. I_num is the numerical value of integration, h is the step size, and E is the error, which is the absolute difference between true and numerical results. First, we are going to determine the true value of the error function while integrating from 0 to 1.5. Then, we will apply the corresponding methods and compare the absolute error.

MATLAB Code for Midpoint Rectangle Rule

close all;
clear;
clc;

% Defining the error function
f = @(x) (2/sqrt(pi))*exp(-x.^2);
a = 0;
b = 1.5; % Upper integration limit
n = [1 2 4 8 16 32 64 128];
L = length(n);

% Exact integral using MATLAB built-in function
I = integral(f,a,b);

fprintf('\n     n          h         I_NUM          E')

for j = 1:L
h(j)=(b-a)/n(j);

x_1 = a:h(j):b;
% Using MATLAB built-in command for repeated arrays
x_2 = repelem(x_1(2:end-1),2);
xh = [a x_2 b];

% Numerical integration by mid-point rectangular rule
I_NUM=0;
for i = 1:length(x_1)-1
I_NUM = I_NUM + f(a  +(i-1/2).*h).*h;
end

% Defining the error
E(j) = abs(I-I_NUM(j));

fprintf('\n %5.6g %11.6g %12.6g %10.6g', n(j), h(j), I_NUM(j), E(j))
fprintf('\n')

end

The following table shows the parameters determined by this approach:

MATLAB Code for Trapezoidal Rule

close all;
clear;
clc;

% Defining the error function
f = @(x) (2/sqrt(pi))*exp(-x.^2);
a = 0;
b = 1.5; % Upper integration limit
n = [1 2 4 8 16 32 64 128];
L = length(n);

% Exact integral using MATLAB built-in function
I = integral(f,a,b);

fprintf('\n     n          h         I_NUM          E')

for j = 1:L
h(j)=(b-a)/n(j);

% Numerical integration using Trapezoidal rule
I_NUM=0;
for i = 1:n(j)
I_NUM = I_NUM+(f(a+(i-1).*h) + f(a+i.*h)).*h/2;
end

% Defining the error
E(j) = I - I_NUM(j);

fprintf('\n %5.6g %11.6g %12.6g %10.6g', n(j), h(j), I_NUM(j), E(j))
fprintf('\n')
end

The following table shows the parameters determined by this approach:

MATLAB Code for Simpson’s 1/3 Rule

close all;
clear;
clc;

% Defining the error function
f = @(x) (2/sqrt(pi))*exp(-x.^2);
a = 0;
b = 1.5; % Upper integration limit
n = [1 2 4 8 16 32 64 128];
L = length(n);

% Exact integral using MATLAB built-in function
I = integral(f,a,b);

fprintf('\n     n          h         I_NUM          E')

for j = 1:L
h(j) = (b-a)/(2*n(j));

% Numerical integration using Simpson’s 1/3 rule

I_NUM = 0;

for i = 1:n(j)
I_NUM = I_NUM+(f(a+(2*i-2).*h) + 4*f(a+(2*i-1).*h) + f(a+(2*i).*h)).*h./3;
end
E(j) = abs(I - I_NUM(j));

fprintf('\n %5.6g %11.6g %12.6g %10.6g', n(j), h(j), I_NUM(j), E(j))
fprintf('\n')
end

The following table shows the parameters determined by this approach:

MATLAB Code for Simpson’s 3/8 Rule

close all;
clear;
clc;

% Defining the error function
f = @(x) (2/sqrt(pi))*exp(-x.^2);
a = 0;
b = 1.5; % Upper integration limit
n = [1 2 4 8 16 32 64 128];
L = length(n);

% Exact integral using MATLAB built-in function
I = integral(f,a,b);

fprintf('\n     n          h         I_NUM          E')

for j = 1:L
h(j) = (b-a)/(3*n(j));

% Numerical integration using Simpson’s 3/8 rule
I_NUM = 0;
for i=1:n(j)
I_NUM = I_NUM + (f(a+(3*i-3).*h) + 3*f(a+(3*i-2).*h) + 3*f(a+(3*i-1).*h) + f(a+(3*i).*h))*3.*h/8;
end

% Defining the error
E(j) = abs(I - I_NUM(j));

fprintf('\n %5.6g %11.6g %12.6g %10.6g', n(j), h(j), I_NUM(j), E(j))
fprintf('\n')

end

The following table shows the parameters determined by this approach:

We see from the above comparisons that both Simpson's rules give us a good approximation of the true value, especially when n is greater than 4. For the rectangular approach, the value reaches close to the true value when n is 128. The same is also evident for the trapezoidal rule. So, we can conclude that both the Simpson's rules are more efficient and accurate approaches compared to the rest.

#TrapezoidalRule #SimpsonsRule #Matlab #NumericalIntegration #Blog #Blogger