Trapezoidal Rule

Numerical Methods

The trapezoidal rule is a simple and widely-used numerical integration method. For periodic functions, it converges exponentially, making it ideal for Fourier coefficient calculations.

Derivation

We start with the integral of a function $f (x)$ over an interval $[a, b]$ :

I = \int_{a}^{b} f (x) d x

Using the Fundamental Theorem of Calculus, we approximate $f (x)$ by its linear interpolation between $x = a$ and $x = b$ :

f (x) \approx f (a) + \frac{x - a}{b - a} (f (b) - f (a))

Substituting this approximation into the integral gives:

I \approx \int_{a}^{b} (f (a) + \frac{x - a}{b - a} (f (b) - f (a))) d x

Splitting the integral into two parts and solving:

I \approx \int_{a}^{b} f (a) d x + \int_{a}^{b} \frac{x - a}{b - a} (f (b) - f (a)) d x

Evaluating the first integral:

\int_{a}^{b} f (a) d x = f (a) (b - a)

For the second integral, we compute:

\int_{a}^{b} \frac{x - a}{b - a} (f (b) - f (a)) d x = (f (b) - f (a)) \frac{( b - a ) ^{2}}{2 ( b - a )} = \frac{( b - a )}{2} (f (b) - f (a))

Summing both parts:

I \approx (b - a) f (a) + \frac{( b - a )}{2} (f (b) - f (a))

Factoring $(b - a)$ :

I \approx \frac{( b - a )}{2} (f (a) + f (b))

This results in the basic trapezoidal rule for a single interval:

\int_{a}^{b} f (x) d x \approx \frac{( b - a )}{2} (f (a) + f (b))

Formula

For an integral:

I = \int_{a}^{b} f (x) d x

The trapezoidal rule with $N$ subintervals is:

I \approx \frac{h}{2} [f (a) + 2 i = 1 \sum N - 1 f (a + ih) + f (b)]

where $h = (b - a) / N$ is the step size.

Error Analysis

The error for the trapezoidal rule is:

E = - \frac{( b - a ) ^{3}}{12 N ^{2}} f^{''} (ξ)

for some $ξ \in [a, b]$ . This gives $O (h^{2})$ convergence for general functions.

Special case: For periodic functions, the trapezoidal rule converges exponentially, with error $O (e^{- c N})$ for some constant $c$ .

C++ Implementation

The following code implements the trapezoidal rule with a template for flexibility:

#include <functional>
#include <cmath>

// Template-based trapezoidal rule
template<typename F>
double trap(F f, double a, double b, int N) {
    double h = (b - a) / N;
    double s = 0.5 * (f(a) + f(b));
    for (int i = 1; i < N; ++i)
        s += f(a + i*h);
    return s * h;
}

// Function-based version
double trapezoidalRule(double a, double b, double dx, 
                       std::function<double(double)> func) {
    size_t points = std::abs(std::ceil((b - a)/dx));
    double result = 0.0;
    
    for (int i = 0; i < points - 1; i++) {
        result += func(a + i*dx) + func(a + (i+1)*dx);
    }
    result = result * dx * 0.5;
    return result;
}

// Example: Integrate Gaussian
void testGaussianIntegration() {
    double alpha = 1.0;
    auto gauss = [alpha](double x) { 
        return std::exp(-alpha * x * x); 
    };
    
    double analytical = std::sqrt(std::numbers::pi_v<double> / alpha);
    
    // Using template version
    double integral1 = trap(gauss, -5.0, 5.0, 10000);
    
    // Using function version
    std::function<double(double)> gauss_func = gauss;
    double integral2 = trapezoidalRule(-5.0, 5.0, 0.001, gauss_func);
    
    std::cout << "Analytical: " << analytical << std::endl;
    std::cout << "Template version: " << integral1 << std::endl;
    std::cout << "Function version: " << integral2 << std::endl;
}

// Example: Periodic function (exponential convergence)
void testPeriodicFunction() {
    auto periodic = [](double x) {
        return std::sin(2.0 * std::numbers::pi_v<double> * x);
    };
    
    // Integrate over one period [0, 1]
    // Exact result is 0
    double integral = trap(periodic, 0.0, 1.0, 100);
    std::cout << "Integral of sin(2πx) over [0,1]: " << integral << std::endl;
    std::cout << "Error: " << std::abs(integral) << std::endl;
    // Error should be very small due to exponential convergence
}

int main() {
    testGaussianIntegration();
    testPeriodicFunction();
    return 0;
}

Key Features

Template version for zero-overhead function calls
Function-based version for runtime flexibility
Exponential convergence for periodic functions
Simple and robust implementation

Comparison with Other Methods

Method	Convergence (General)	Convergence (Periodic)
Trapezoidal	$O (h^{2})$	$O (e^{- c N})$
Simpson's	$O (h^{4})$	$O (e^{- c N})$
Monte Carlo	$O (1/ N)$	$O (1/ N)$

Worked examples

Estimate an integral with the trapezoidal rule → Approximate ∫₀¹ e^(−x²) dx using the trapezoidal rule with 4 sub-intervals, then bound the error and compare to the true value.

Browse all worked examples →