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

### 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