Computational Physics

A6 Numerical quadratures 567

Such an approximation is useless without an estimate of the error in the result. Iffis
continuous and bounded on the interval[xn,xn+h], its value does not deviate from
f(xn)more thanMh, withMsome finite constant, depending on the functionfand
the values ofaandb. Therefore, the relative error in this result isO(h): doubling
the number of integration points reduces the error by a factor of about 2.
A way of approximating the integral of a differentiable functionfto second order
inhusing the grid points given in Eq. (A.25) is the following. Approximatef on
the interval[xn,xn+ 1 ]to first order inhby

f(xn)+

(x−xn)

h

[f(xn+ 1 )−f(xn)]+O(h^2 ),n=0,...,N−1. (A.27)

Integrating this form analytically on the interval, we obtain
∫b

a

f(x)dx=h

[

1

2

f(x 0 )+f(x 1 )+f(x 2 )+···+f(xN− 1 )+

1

2

f(xN)

]

+O(h^2 ).
(A.28)
This is called the trapezoidal rule.
Iff is twice differentiable, we can approximate it by piecewise quadratic
functions. For a single interval using three equidistant points, we have
∫b

a

f(x)dx=

b−a

6

[

f(a)+ 4 f

(

a+b

2

)

+f(b)

]

+O[(b−a)^5 ], (A.29)

which is one order of magnitude better than expected, because there is a cancellation
of errors at the left and right hand boundaries, resulting from the symmetry of this
formula. If the interval[a,b]is split up intoN/2 subintervals, each of size 2h,we
obtain the following expression:
∫b

a

f(x)dx=h

1

3

[f(x 0 )+ 4 f(x 1 )+ 2 f(x 3 )+ 4 f(x 4 )+···

+ 2 f(xN− 2 )+ 4 f(xN− 1 )+f(xN)]+O(h^4 ). (A.30)

This is an example ofSimpson’s rule.
It would be easy to extend this list of algorithms. They all boil down to making
some polynomial approximation of the integrand and then integrating the approx-
imate formula analytically, which is trivial for a piecewise polynomial. Also, many
tricks exist for integrating functions containing singularities, but these are beyond
the scope of this discussion. In all cases, the order of the accuracy is equal to the
number of points used in fitting the polynomial, since fornpoints one can determine
a polynomial of ordern−1 having the same value as the function to be integrated
on allnpoints.
A very efficient method consists of repeating the trapezoidal rule and performing
it for successive values ofh, each having half the size of the previous one. This
yields a sequence of approximants to the integral for various values ofh. These