Numerical integration is the process of approximating a definite integral when an analytical solution is difficult to obtain or when the function is available only through discrete samples.

Consider

\[I = \int_{t_0}^{t_f} L(t) \, dt.\]

Numerical integration is widely used in engineering and scientific applications. For example, integrating velocity gives displacement, integrating power gives energy, and integrating acceleration gives the change in velocity. It is also useful when processing experimental data, where a quantity is known only at discrete time instants.

In many practical problems, we therefore have a set of values

\[L_0, L_1, \dots, L_N\]

rather than an analytical expression for $L(t)$. The basic idea is to approximate the function between these points using a simple polynomial and then integrate that polynomial.

Two commonly used methods are the trapezoidal rule, which uses a linear approximation, and Simpson’s rule, which uses a quadratic approximation.

Time domain

Consider the interval

\[t \in [t_0, t_f].\]

Divide it into $N$ equally spaced intervals:

\[t_0, t_1, \dots, t_N.\]

The spacing between consecutive points is

\[h = t_{i+1} - t_i = \frac{t_f - t_0}{N}.\]

Therefore,

\[t_i = t_0 + i h.\]

Let

\[L_i = L(t_i).\]

The problem is now to approximate

\[I = \int_{t_0}^{t_f} L(t) \, dt\]

using the discrete values

\[L_0, L_1, \dots, L_N.\]

The two methods differ in how they approximate $L(t)$ between the sampled points.

Trapezoidal method

Basic idea

Consider one interval $[t_i, t_{i+1}]$. Approximate $L(t)$ by the straight line joining

\[(t_i, L_i) \quad \text{and} \quad (t_{i+1}, L_{i+1}).\]

The area under this linear approximation is the area of a trapezoid:

\[\int_{t_i}^{t_{i+1}} L(t) \, dt \approx \frac{h}{2} (L_i + L_{i+1}).\]

This is the trapezoidal rule.

Derivation

Introduce a local coordinate

\[\tau = t - t_i, \quad 0 \le \tau \le h.\]

The linear approximation is

\[L_i(\tau) = L_i + \frac{L_{i+1} - L_i}{h} \tau.\]

Integrating over the interval,

\[\begin{aligned} \int_{0}^{h} L_i(\tau) \, d\tau &= \int_{0}^{h} \left( L_i + \frac{L_{i+1} - L_i}{h} \tau \right) d\tau \\ &= L_i h + \frac{L_{i+1} - L_i}{h} \frac{h^2}{2} \\ &= \frac{h}{2} (L_i + L_{i+1}). \end{aligned}\]

Thus,

\[\int_{t_i}^{t_{i+1}} L(t) \, dt \approx \frac{h}{2} (L_i + L_{i+1}).\]

Composite trapezoidal rule

The same approximation can be applied to every interval:

\[\int_{t_0}^{t_f} L(t) \, dt \approx \sum_{i=0}^{N-1} \frac{h}{2} (L_i + L_{i+1}).\]

Collecting the terms gives

\[\int_{t_0}^{t_f} L(t) \, dt \approx h \left[ \frac{L_0 + L_N}{2} + \sum_{i=1}^{N-1} L_i \right].\]

Thus, the composite trapezoidal rule approximates the curve by a sequence of straight-line segments.

Simpson method

Basic idea

The trapezoidal rule uses a linear approximation. We can obtain a higher-order approximation by fitting a quadratic polynomial through three points.

Consider two consecutive intervals,

\[t_i, t_{i+1}, t_{i+2},\]

with equal spacing

\[t_{i+1} - t_i = t_{i+2} - t_{i+1} = h.\]

A quadratic polynomial is fitted through these three points.

Derivation

Introduce the local coordinate

\[\tau = t - t_i, \quad 0 \le \tau \le 2h.\]

Approximate $L(t)$ by

\[L_i(\tau) = a \tau^2 + b \tau + c.\]

The polynomial must pass through the three sampled values:

\[L_i(0) = L_i, \quad L_i(h) = L_{i+1}, \quad \text{and} \quad L_i(2h) = L_{i+2}.\]

Therefore,

\[c = L_i,\] \[a h^2 + b h = L_{i+1} - L_i,\] \[4 a h^2 + 2 b h = L_{i+2} - L_i.\]

Solving these equations,

\[a = \frac{L_{i+2} - 2 L_{i+1} + L_i}{2 h^2},\] \[b = \frac{-L_{i+2} + 4 L_{i+1} - 3 L_i}{2 h}.\]

The integral of the quadratic approximation is

\[\int_{0}^{2h} L_i(\tau) \, d\tau = \frac{8 a h^3}{3} + 2 b h^2 + 2 c h.\]

Substituting $a$, $b$, and $c = L_i$ gives

\[\int_{t_i}^{t_{i+2}} L(t) \, dt \approx \frac{h}{3} (L_i + 4 L_{i+1} + L_{i+2}).\]

This is Simpson’s rule.

Composite Simpson rule

Simpson’s rule is applied over pairs of intervals. Therefore, $N$ must be even.

Summing the individual panels gives

\[\int_{t_0}^{t_f} L(t) \, dt \approx \frac{h}{3} \left[ L_0 + L_N + 4 \sum_{\substack{i=1 \\ i \text{ odd}}}^{N-1} L_i + 2 \sum_{\substack{i=2 \\ i \text{ even}}}^{N-2} L_i \right].\]

The resulting weights are

\[1, \; 4, \; 2, \; 4, \; 2, \; \dots, \; 4, \; 1.\]

For example, for $N=6$,

\[\int_{t_0}^{t_6} L(t) \, dt \approx \frac{h}{3} \left[ L_0 + 4 L_1 + 2 L_2 + 4 L_3 + 2 L_4 + 4 L_5 + L_6 \right].\]

Comparison

The two methods can be summarized as follows:

Method Approximation Intervals per panel Global error
Trapezoidal Linear 1 $\mathcal{O}(h^2)$
Simpson Quadratic 2 $\mathcal{O}(h^4)$

The trapezoidal rule connects neighboring samples using straight lines, while Simpson’s rule fits a parabola through every three consecutive samples.

For a sufficiently smooth function, the composite trapezoidal rule has a global error of order $\mathcal{O}(h^2)$, while the composite Simpson rule has a global error of order $\mathcal{O}(h^4)$.

Simpson’s rule therefore provides a higher-order approximation, but it requires an even number of intervals.

Limitations

It may seem that increasing the polynomial degree should always improve the approximation. However, using a single high-degree polynomial over many uniformly spaced points can lead to the Runge phenomenon and numerical ill-conditioning.

For this reason, numerical integration generally uses piecewise low-degree polynomial approximations rather than a single high-degree polynomial over the entire interval.

Higher-order integration methods can instead use carefully chosen non-uniform points or other approximation techniques, such as Gaussian quadrature.

Numerical code

The trapezoidal and Simpson quadrature is demonstrated for the integral $ \int_{0}^{2} \sin (x)=2$ in python. The example requires the numpy package for the math and matplotlib for visualization of error plots.

import numpy as np
import matplotlib.pyplot as plt


def L(t):
    return np.sin(t)


exact = 2.0
t0 = 0.0
tf = np.pi


def trapezoidal(f, t0, tf, N):
    h = (tf - t0) / N
    t = np.linspace(t0, tf, N + 1)
    y = f(t)
    return h * (0.5 * y[0] + np.sum(y[1:-1]) + 0.5 * y[-1])


def simpson(f, t0, tf, N):
    if N % 2 != 0:
        raise ValueError("N must be even for Simpson's rule.")
    h = (tf - t0) / N
    t = np.linspace(t0, tf, N + 1)
    y = f(t)
    return (h / 3) * (y[0] + y[-1] + 4 * np.sum(y[1:-1:2]) + 2 * np.sum(y[2:-2:2]))


N_values = np.array([4, 8, 16, 32, 64, 128, 256, 512])
h_values = (tf - t0) / N_values

trap_values = [trapezoidal(L, t0, tf, N) for N in N_values]
simp_values = [simpson(L, t0, tf, N) for N in N_values]

trap_errors = [abs(I_trap - exact) for I_trap in trap_values]
simp_errors = [abs(I_simp - exact) for I_simp in simp_values]

plt.figure(figsize=(7, 5))

plt.loglog(h_values, trap_errors, "o-", label="Trapezoidal")
plt.loglog(h_values, simp_errors, "s-", label="Simpson")

plt.xlabel(r"$h$")
plt.ylabel("Absolute error")
plt.title(r"Numerical integration of $\sin(t)$ over $[0,\pi]$")
plt.grid(True, which="both", linestyle="--", alpha=0.5)
plt.legend()

plt.tight_layout()
plt.show()

Conclusion

In this post, we derived and implemented two fundamental numerical integration techniques—the trapezoidal and Simpson’s rules in Python. By testing them on $\sin(t)$ over $[0, \pi]$, we verified that Simpson’s quadratic approximation achieves significantly higher accuracy than the trapezoidal rule for the same step size, and we visualized their second- and fourth-order convergence rates on a log-log error plot.