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

## No comments:

## Post a Comment