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