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