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.

First Order Ordinary Differential Equations









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.

Defining differential equations in MATLAB












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.

Application of the explicit Euler method in MATLAB
















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.

Final stage where we call the function



 











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.

Plot showing the solutions of differential equations














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).

Error function integral




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:

Table showing integration 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:
Show comparison with true integration value












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:
Show comparison with true integration value












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:
Showing comparison with true integration value












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:
Showing comparison with true integration value













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