MATLAB - Trapezoidal Rule

MATLAB - Numerical Integration with the Trapezoidal Rule

Numerical integration is a fundamental technique for approximating the definite integral of a function, especially when an analytical solution is difficult or impossible to find. The trapezoidal rule is a popular numerical method that approximates the area under a curve by dividing it into a series of trapezoids and summing their areas. This approach is often more accurate than using rectangles.

The Theory Behind the Trapezoidal Rule

Step 1: Divide the Interval [a, b]

First, we divide the integration interval [a, b] into n smaller subintervals, each of equal width, Δx.

The width Δx is calculated as: $$ \Delta x = \frac{b - a}{n} $$

This creates n+1 points: x₀, x₁, x₂, …, xₙ, where:

$$ x_i = a + i \cdot \Delta x \quad \text{for } i = 0, 1, 2, \dots, n $$

Here, x₀ = a and xₙ = b are the endpoints of the interval.

Step 2: Form Trapezoids and Calculate Their Areas

For each subinterval [xᵢ₋₁, xᵢ], we form a trapezoid. The parallel sides of the trapezoid are the function values at the endpoints, f(xᵢ₋₁) and f(xᵢ), and its height is Δx.

The area of a single trapezoid is given by the standard formula:

$$ \text{Area}i = \frac{1}{2} \times (\text{sum of parallel sides}) \times \text{height} = \frac{\Delta x}{2} [f(x{i-1}) + f(x_{i})] $$

Step 3: Sum the Areas

The total approximate integral is the sum of the areas of all n trapezoids:

$$ \int_{a}^{b} f(x) , dx \approx \sum_{i=1}^{n} \text{Area}_i $$

Expanding this sum reveals a more efficient formula. Notice that all interior points (f(x₁) through f(xₙ₋₁)) are shared by two adjacent trapezoids, so they are counted twice.

$$ \int_{a}^{b} f(x) , dx \approx \frac{\Delta x}{2} \left[ f(x_0) + 2f(x_1) + 2f(x_2) + \dots + 2f(x_{n-1}) + f(x_n) \right] $$

Implementing the Trapezoidal Rule in MATLAB

MATLAB provides the built-in trapz function to perform trapezoidal integration directly on discrete data points. This is extremely useful when you have experimental data (e.g., sensor readings over time) rather than a mathematical function.

Syntax

Q = trapz(Y)
Q = trapz(X, Y)
Q = trapz(___, dim)

Syntax Explanation

Q = trapz(Y): Computes the integral of the data in vector Y assuming unit spacing (i.e., X = [1, 2, 3, ...]). If Y is a matrix, trapz integrates each column independently.

Q = trapz(X, Y): Computes the integral of Y with respect to the coordinates specified in X. X and Y must be vectors of the same length. This is the most common usage for non-unit spacing.

Q = trapz(___, dim): Integrates along a specified dimension dim for multi-dimensional arrays. For example, trapz(X, Y, 2) integrates along the rows of matrix Y.

Example 1: Integrating Data with Uniform and Non-Uniform Spacing

Let’s calculate the integral of the function y = x² from x=1 to x=5 using discrete points.

% Define the coordinate and value vectors
X = 1:5; % X = [1, 2, 3, 4, 5]
Y = X.^2; % Y = [1, 4, 9, 16, 25]

% Case 1: Using trapz(Y) assumes unit spacing (which is correct here)
Q1 = trapz(Y);

% Case 2: Using trapz(X, Y) explicitly provides the spacing
Q2 = trapz(X, Y);

% The analytical solution is ∫x² dx from 1 to 5 = [x³/3] from 1 to 5 = 125/3 - 1/3 = 124/3 ≈ 41.33
fprintf('Result with assumed unit spacing (Q1): %.4f\n', Q1);
fprintf('Result with explicit spacing (Q2): %.4f\n', Q2);
fprintf('Analytical Result: %.4f\n', 124/3);

% --- Output ---
% Result with assumed unit spacing (Q1): 42.0000
% Result with explicit spacing (Q2): 42.0000
% Analytical Result: 41.3333

The numerical result 42.0 is a good approximation of the true value 41.333. The accuracy can be improved by using more data points.

Example 2: Integrating a Matrix Along Columns and Rows

% Create a 3x3 matrix representing 3 different datasets
Y = [1 2 3;
     4 5 6;
     7 8 9];

% Integrate along each COLUMN (dim=1 is the default)
Q_cols = trapz(Y);

% Integrate along each ROW (specify dim=2)
Q_rows = trapz(Y, 2);

disp('Original Matrix Y:');
disp(Y);
disp('Integral of each column:');
disp(Q_cols);
disp('Integral of each row:');
disp(Q_rows);

% --- Output ---
% Original Matrix Y:
%      1     2     3
%      4     5     6
%      7     8     9
%
% Integral of each column:
%      8    10    12
%
% Integral of each row:
%      4
%     10
%     16

Practical Application: Calculating Distance from Velocity Data

Imagine a sensor records a vehicle’s velocity at different time points. We can integrate the velocity to find the total distance traveled.

% Time in seconds (non-uniform spacing)
time_s = [0, 10, 25, 40, 60];

% Velocity in meters/second
velocity_mps = [0, 22, 28, 30, 25];

% Calculate total distance by integrating velocity over time
total_distance_m = trapz(time_s, velocity_mps);

fprintf('Total distance traveled: %.2f meters\n', total_distance_m);

% --- Output ---
% Total distance traveled: 1475.00 meters

Modern Alternative: The integral Function

For situations where you have a mathematical function to integrate, MATLAB’s integral function is the modern, more accurate, and preferred choice. It uses a more sophisticated technique called adaptive quadrature.

% Define the function to integrate using a function handle
fun = @(x) x.^2;

% Define the integration interval
a = 1;
b = 5;

% Calculate the integral
Q = integral(fun, a, b);

fprintf('The result from integral() is: %.4f\n', Q);

% --- Output ---
% The result from integral() is: 41.3333

Best Practices and Common Pitfalls

• trapz vs. integral: Use trapz when you have discrete data points (e.g., from an experiment). Use integral when you have a symbolic function to evaluate.
• Monotonic Coordinates: When using trapz(X,Y), ensure the coordinate vector X is monotonic (either always increasing or always decreasing).
• Vector Length: X and Y must have the same number of elements.
• Accuracy: The accuracy of the trapezoidal rule depends on the number of data points. More points (a smaller Δx) generally lead to a more accurate result.