Trapezoidal method, also known as trapezium method or simply trapezoidal rule, is a popular method for numerical integration of various functions (approximation of definite integrals) that arise in science and engineering. This method is mainly applicable to estimate the area under a curve by splitting the entire area into a number of trapeziums of known area.
In this tutorial, we are going to write a program code for Trapezoidal method in MATLAB, going through its mathematical derivation and a numerical example. You can check out our earlier tutorials where we discussed a C program and algorithm/flowchart for this method.
Derivation of Trapezoidal Method:
Consider a curve f(x) = 0 as shown below:
Here, the area under the curve f(x) from x0 to xn is to be determined which is theoretically nothing but the numerical integration of curve f(x) in the interval [x0, xn]. So, in order to estimate the area, we split the whole area into ‘n’ number of trapeziums or trapezoids of equal height ‘h’.
Let y0, y1, y2….. yn be the ordinates of the curve at various abscissas.
The area of first trapezium, a1 = h * (y0 +y1)/ 2
The area of second trapezium, a2 = h * (y1 +y2)
The area of third trapezium, a3 = h * (y2 +y3)/ 2
… . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . .. . .. . .
. . . .. . . . .. . . . .. . . . .. . . . . . .. .. . . . .. . . . .. . .
The area of nth trapezium, an = h * (yn-1 +yn)/ 2
The total area under the curve is the summation of all these small areas.
So, A = a1 + a2 + a3 + ….. + an
or, A= h * (y0 +y1)/ 2 + h * (y1 +y2) + h * (y2 +y3)/ 2+ ….. + h * (yn-1 +yn)/ 2
or, A= h { (y0 + yn)/2 + ( y1 + y2+ …….. + yn-1)}
Note: This formula is the required expression which will be used in the program code for Trapezoidal rule in MATLAB.
Alternatively, it can be concluded that:
Trapezoidal Method in MATLAB:
>> n = 5; % number of small trapeziums formed after splitting
a = 1.0; % starting point or lower limit of the area
b = 2.0; % end point or upper limit of the area
sum = 0.0; % to find the sum
dx = (b-a)/(n-1); % to find step size or height of trapezium
% Generating the samples
for i = 1:n
x(i) = a + (i-1)*dx;
end
% to generate the value of function at different values of x or sample
for i = 1:n
y(i) = x(i).^2;
end
% computation of area by using the technique of trapezium method
for i = 1:n
if ( i == 1 || i == n) % for finding the sum of fist and last ordinate
sum = sum + y(i)./2;
else
sum = sum + y(i); % for calculating the sum of other ordinates
end
end
area = sum * dx; % defining the area
area
The above code for Trapezoidal method in MATLAB has been programmed to find the area under the curve f(x) = x2 in the interval [1, 2]. So, the user doesn’t need to give any input to the program.
In this program, the whole area has been divided into five trapeziums. The program calculates the area of each trapezium as explained in the derivation and finds the required sum.
If it is intended to find the area under other curves using this code, the user has to change the limits and the value of the function in the source code.
To run the program, copy the given code in MATLAB command window and hit enter. Here’s a screenshot sample of the output you’ll get:
As the aforementioned MATLAB code has an inbuilt function, it may be difficult (for some of you) to modify the source code for other functions and limits. In such case, you can use this program code instead:
function I = trapezoidal_f1 ( f )
% for calculating integrals using trapezoidal rule when function is known
%asking for the range and desired accuracy
R= input( ' Enter the limits of integrations [ x_min, x_max] :\n');
tol = input(' Error allowed in the final answer should be of an order : \n');
a= R(1,1); b = R(1,2);
%initial h and n
n = 100;
h = (b -a )/100;
%for calculating maximum of h^2 *f''(x) in the given region
for k = 0:100
x( 1, k+1 ) = a + k*h ;
y2 ( 1, k+1) = feval ( f, x(1,k+1) + 2*h ) - 2*feval( f, x(1,k+1) + h ) + feval( f, x(k+1) );
end
[ y i ] = max( y2);
x_opt = x(1,i);
% for calculating the desired value of h
m=0;
while abs((feval ( f, x_opt + 2*h ) - 2*feval(f, x_opt+ h) + feval(f, x_opt)) * ( b-a)/12) > tol % Global error for trapezoidal rule
m = m +1;
h = h * 10^-m;
end
%calculating n
n = ceil( (b-a)/h );
h = ( b-a )/ n;
% forming matrix X
for k = 1:(n+1)
X(k,1) = a + (k-1)*h;
X(k,2) = feval( f, X(k,1));
end
% trapezoidal formula
I = h/2 * (2*sum( X(:,2))- X(1,2)- X( n,2));
% displaying final result
h
n
disp(sprintf(' The area under this curve in the interval ( %10f , %10f ) is %10.6f .',a,b,I))
When this MATLAB program is executed, it asks for the value of the interval in which the area of curve is to be estimated. After that, the user needs to input allowed error in calculation. Finally, the area under the curve is displayed as output.
How to input the function?
To run the code, follow these steps:
- Copy the source code given above in MATLAB editor and save as file_mam file.
- Create another MATLAB file in the same folder in which .m is saved.
- In this newly created MATLAB file, define the function following the syntax given below.
function I = trapezoidal_f1 ( f )
f = x^2
The sample output of this program is given below:
Trapezoidal Method Numerical Example:
Now, let’s analyze numerically the first program for Trapezoidal method in MATLAB. We have to find the area under the curve f(x) = x2 in the interval [1, 2].
Given function, f(x) = x2
x0 = 1 and xn = 2
Take number of steps, n = 5 such that the height of the interval, h = ( 2-1)/5 = 0.2. Detailed calculation is presented in the table below:
Now, using the formula for trapezoidal method:
Area = 0.2/2 * [ (1+4) + 2*( 1.44 + 1.96 + 2.56 + 3.24)] = 2.3438
Thus, the answer is same as the one obtained using the program for Trapezoidal method in MATLAB.
If you have any questions regarding trapezoidal method or its MATLAB code, bring them up to me from the comments box below. Also, you can find more Numerical Methods tutorials here.
