MATLAB - Simpson's Rule

MATLAB - Numerical Integration with Simpson’s Rule

Simpson’s Rule is a powerful numerical method used to approximate the value of a definite integral. It offers greater accuracy than simpler methods like the Trapezoidal Rule by using quadratic polynomials to approximate the function within each subinterval. This makes it particularly effective for functions that are smooth and well-behaved.

The rule is most commonly applied as Simpson’s 1/3 Rule, which requires the integration interval [a, b] to be divided into an even number (n) of subintervals.

The formula for Simpson’s 1/3 Rule is:

$$\int_{a}^{b} f(x) , dx \approx \frac{\Delta x}{3} [f(x_{0}) + 4f(x_{1}) + 2f(x_{2}) + \dots + 4f(x_{n-1}) + f(x_{n})]$$

where $$\Delta x = \frac{b-a}{n}$$ and $$x_i = a + i\Delta x$$. Notice the pattern of weights: 1, 4, 2, 4, ..., 2, 4, 1.

How Simpson’s Rule Works

1. Divide the Interval: The total interval [a, b] is split into n equal subintervals, where n must be even.
2. Fit Quadratic Polynomials: Across each pair of adjacent subintervals, a unique quadratic polynomial (a parabola) is fitted to pass through the three points of the function.
3. Sum the Areas: The integral is approximated by summing the exact areas under each of these parabolas.

Implementation in MATLAB

While you can implement Simpson’s Rule from scratch to understand the algorithm, modern MATLAB offers more efficient methods. We will explore three approaches: a basic loop, a vectorized implementation, and using MATLAB’s built-in function.

Approach 1: Implementation with a for Loop (for learning)

This method is straightforward and clearly shows the algorithm’s logic. We will approximate the integral of $$f(x) = x^2$$ from a=0 to b=2.

% 1. Define the function, interval, and number of subintervals
f = @(x) x.^2;
a = 0;
b = 2;
n = 10; % Must be an even number

% 2. Check if n is even
if mod(n, 2) ~= 0
    error('Number of subintervals (n) must be even for Simpson''s 1/3 rule.');
end

% 3. Calculate step size
h = (b - a) / n;

% 4. Initialize sum with the first and last terms
integral_sum = f(a) + f(b);

% 5. Loop through interior points to add weighted terms
for i = 1:(n - 1)
    x_i = a + i * h;
    if mod(i, 2) == 0 % Even indices get weight 2
        integral_sum = integral_sum + 2 * f(x_i);
    else % Odd indices get weight 4
        integral_sum = integral_sum + 4 * f(x_i);
    end
end

% 6. Finalize the calculation
approx_integral = (h / 3) * integral_sum;

fprintf('Approximate integral (loop): %.6f\n', approx_integral);

% For comparison, the exact value is 8/3 ≈ 2.666667

Output: Approximate integral (loop): 2.666667

Approach 2: Vectorized Implementation (Best Practice)

A more efficient and ‘MATLAB-style’ approach avoids loops by using vector operations. This is significantly faster for large n.

% 1. Define function, interval, and n (as before)
f = @(x) x.^2; a = 0; b = 2; n = 10;
if mod(n, 2) ~= 0, error('n must be even'); end

% 2. Create a vector of all x points
x = linspace(a, b, n + 1);

% 3. Calculate the corresponding y values
y = f(x);

% 4. Apply Simpson's rule using vector summation
% Sum odd-indexed interior points (y(2), y(4), ...)
% Sum even-indexed interior points (y(3), y(5), ...)
h = (b - a) / n;
approx_integral_vec = (h/3) * (y(1) + 4*sum(y(2:2:n)) + 2*sum(y(3:2:n-1)) + y(n+1));

fprintf('Approximate integral (vectorized): %.6f\n', approx_integral_vec);

Output: Approximate integral (vectorized): 2.666667

Approach 3: The Modern MATLAB Way with integral (Recommended)

For most practical applications, you should use MATLAB’s built-in integral function. It uses a globally adaptive quadrature, which is more robust and accurate than a fixed-step rule like Simpson’s. It automatically handles complexities and provides an error estimate.

% 1. Define the function and interval (as before)
f = @(x) x.^2;
a = 0;
b = 2;

% 2. Call the integral function
matlab_integral = integral(f, a, b);

fprintf('Integral using built-in function: %.6f\n', matlab_integral);

Output: Integral using built-in function: 2.666667

Conclusion: When to Use Each Method

• Loop Implementation: Use for learning the algorithm or when you need to inspect intermediate steps for educational purposes.
• Vectorized Implementation: Use when you need to implement Simpson’s rule specifically (e.g., for a class assignment) and want an efficient, well-written script.
• integral Function: Use for all practical and professional work. It is highly optimized, more accurate, and can handle a wider variety of functions, including those with singularities.