Skip to content
Numerical Calculus

Numerical integration: trapezoidal, Simpson, Romberg and Gauss

Trapezoidal and Simpson's rules with their errors, a numerical example, Richardson extrapolation, adaptive integration, Gaussian quadrature and sampled signals.

4 min read 0 comments

Many integrals have no antiderivative expressible with elementary functions (\int e^{-x^2}dx, \int\frac{\sin x}{x}dx, the perimeter of an ellipse), and often the function is known only at points, like the samples of a signal. Numerical integration (or quadrature) approximates the integral by a weighted sum of the function’s values at a few points. It is how a multimeter computes the RMS value of a voltage, how an electronic meter computes the energy used and how the area under an experimental curve is estimated.

The common idea

Split [a, b] into n intervals of width h = (b - a)/n, with nodes x_i = a + ih, approximate f on each interval with a polynomial and integrate the polynomial exactly.

Numerical integration with the trapezoidal rule and with Simpson's rule
Left, the trapezoidal rule: the function is replaced by segments between the nodes. Right, Simpson’s rule: on each pair of intervals the function is replaced by a parabola through three points, which follows it much better.

The trapezoidal rule

On each interval f is replaced by the segment between the two ends, and the area is that of a trapezium. Adding up:

    \[ \int_a^b f\,dx \approx T_h = h\left[\frac{f_0}{2} + f_1 + f_2 + \dots + f_{n-1} + \frac{f_n}{2}\right] \]

The error is

    \[ \int_a^b f\,dx - T_h = -\frac{(b - a)\,h^2}{12}\,f''(\xi) \]

proportional to h^2: halving the step cuts the error by four. The rule is exact for first-degree polynomials. The sign of the second derivative tells which side the error is on: for a convex function the trapezia lie above the curve and overestimate the integral.

Simpson’s rule

With n even, on each pair of intervals f is replaced by the parabola through the three points. Integrating the parabola gives the weights \frac h3(1, 4, 1), and adding up:

    \[ \int_a^b f\,dx \approx S_h = \frac h3\Big[f_0 + 4f_1 + 2f_2 + 4f_3 + \dots + 2f_{n-2} + 4f_{n-1} + f_n\Big] \]

The error is

    \[ -\frac{(b - a)\,h^4}{180}\,f^{(4)}(\xi) \]

proportional to h^4: halving the step cuts the error by sixteen. Surprisingly the rule is exact for cubic polynomials too, not just parabolas: the error of the cubic term cancels by symmetry.

Example. \int_0^1 e^{-x^2}dx = 0.746824\ldots With n = 4 (h = 0.25):

MethodValueError
trapezoidal0.7429843.8\times 10^{-3}
Simpson0.7468553.1\times 10^{-5}

With the same five function values, Simpson is more than a hundred times more accurate.

Proof: Simpson from the trapezoidal rule

The trapezoidal error has the form T_h = I + ch^2 + O(h^4). Computing with steps h and 2h and combining to remove the h^2 term:

    \[ \frac{4T_h - T_{2h}}{3} = I + O(h^4) \]

and working it out gives exactly Simpson’s rule. This is Richardson extrapolation: knowing how the error behaves, you can remove it. Repeating the idea gives Romberg’s method, which reaches high accuracy with few values.

Estimating the error without the derivative

The error formulas contain unknown derivatives. In practice the result is compared at two different steps: for the trapezoidal rule, I - T_h \approx \frac{T_h - T_{2h}}{3}. Adaptive integration programs do this on each subinterval, and add points only where the function changes quickly: few points where it is flat, many near a peak.

Gaussian quadrature

The trapezoidal and Simpson rules use equally spaced nodes. If the node positions can be chosen too, n points can integrate a polynomial of degree 2n - 1 exactly: the Gauss-Legendre formulas. With just two points at \pm\frac{1}{\sqrt3} on [-1, 1]:

    \[ \int_{-1}^{1}f\,dx \approx f\!\left(-\tfrac{1}{\sqrt3}\right) + f\!\left(\tfrac{1}{\sqrt3}\right) \]

exact for every cubic polynomial, like Simpson but with one point fewer. The nodes are the roots of the Legendre polynomials. These formulas are used in finite element programs, where every function evaluation is expensive.

Integrals of sampled signals

For a signal sampled with period T the nodes are fixed by the converter, and the trapezoidal rule or Simpson is applied to the samples. Example: the RMS value of a voltage,

    \[ V_{rms} = \sqrt{\frac{1}{T_p}\int_0^{T_p}v^2(t)\,dt} \]

is computed in “true RMS” multimeters by summing the squared samples over a period. For a periodic signal sampled over a whole number of periods, the trapezoidal rule is surprisingly accurate: the error falls faster than any power of h, because the end errors cancel. It is the same reason the DFT computes the Fourier coefficients of band-limited signals exactly.

Improper integrals and singularities

On infinite intervals the interval is truncated where the function is negligible, or a change of variable makes it finite. If the function has an integrable singularity, such as \frac{1}{\sqrt x} at 0, formulas with nodes at the ends fail (f(0) is infinite) and convergence is slow: it pays to remove the singularity with a change of variable (here x = t^2) before integrating.

Related articles

Subscribe
Notify of
guest
0 Comments
Oldest
Newest Most Voted
0
Would love your thoughts, please comment.x
()
x