Quadrature

Problem: To find an approximate value of $\int _a^bf(x)dx$ when it cannot be evaluated analytically.

 

Figure: Numerical Evaluation of Integrals

The value of such an integral is simply the area under its graph, so we can approximate this by the area of the trapezium with the same corners. This is just $(b-a)\frac{f(a)+f(b)}{2}$ in which we put b-a=h to obtain the trapezium rule $\frac{h}{2}\{ f(a)+f(b)\}$. Clearly there will be a big error if f is strongly curved.

How can we reduce this error? By subdividing the interval of integration into smaller intervals; suppose we divide into 4 equal intervals of width h=(b-a)/4, then the integral is approximated by

$$h\{ \frac{f(a)+f(a+h)}{2}+\frac{f(a+h)+f(a+2h)}{2}+ \frac{f(a+2h)+f(a+3h)}{2}+\frac{f(a+3h)+f(b)}{2}\} .$$

which we rearrange to

$$h\{ \frac{f(a)+f(b)}{2}+f(a+h)+f(a+2h)+f(a+3h)\} .$$

For n intervals h=(b-a)/n and we have

$$\int _a^bf(x)dx\approx h\left\{ \frac{f(a)+f(b)}{2}+\sum _{i=1}^{n-1}f(a+ih)\right\},$$

the composite trapezium rule.

How else might we improve accuracy? By using a curved line to approximate f. We could find a parabola or higher degree polynomial interpolating f at a suitable number of points and integrate this in place of f. We shall achieve the same result, but without evaluating the interpolating polynomial, by the use of the method of undetermined coefficients.

Trapezium Rule
Let $\alpha f(a)+\beta f(b)$ be the quadrature formula for the evaluation of $\int _a^bf(x)dx$. We choose $\alpha ,\beta$ to make the formula exact when $f(x)=1,x,x^2,\ldots$.

$$\begin{array}{lrcl} f(x)=1: & b-a & = & \alpha +\beta \\ f(x)=x: & \frac{1}{2}b^2-a^2 & = & a\alpha +b\beta . \end{array}$$

There is no point in taking further powers of x, since we already have two linear equations in two unknowns. Their solution is $\alpha =\beta =\frac{1}{2} (b-a)$, so the rule is

$$\int _a^bf(x)dx\approx \frac{1}{2}(b-a)\{ f(a)+f(b)\} .$$

Putting a=0,b=h gives the rule in neater form as

$$\int _0^hf(x)dx \approx \frac{1}{2}h\{f(0)+f(h)\} .$$

Simpson's Rule
We wish to find a quadrature rule for $\int _a^bf(x)dx$ which uses the values of f at $a,\frac{a+b}{2},b$. The algebra is much easier if we change the interval [a,b] into an interval which is symmetric about the origin, [-h,h] say. This is achieved by the substitution $x=\frac{b-a}{2}\frac{t}{h}+\frac{a+b}{2}$. Thus

$$\begin{eqnarray*}\int _a^bf(x)dx & = & \frac{b-a}{2h}\int _{-h}^hf(\frac{b-a}{2}... ...{h}+\frac{a+b}{2})dt \\ & = & \frac{b-a}{2h}\int _{-h}^hg(t)dt, \end{eqnarray*}$$

say. We now find a quadrature formula for this integral, which we can insert in order to find a formula for the original integral.

Let

$$\int _{-h}^hg(x)dx\approx \alpha g(-h)+\beta g(0)+\gamma g(h).$$

Using the method of undetermined coefficients, we find

$$\begin{array}{lrcl} g(x)=1: & 2h & = & \alpha +\beta +\gamma ... ...x^2: & \frac{2}{3}h^3 & = & \alpha h^2+\gamma h^2, \end{array}$$

giving $\alpha =\gamma = \frac{h}{3},\beta =\frac{4h}{3}$ and Simpson's rule in the form

$$\int _{-h}^hg(x)dx\approx \frac{h}{3}\{ g(-h)+4g(0)+g(h)\} .$$

Substituting this into the above, we obtain

$$\int _a^bf(x)dx\approx \frac{b-a}{6}\{ f(a)+4f(\frac{a+b}{2})+f(b)\} .$$

NOTE that Simpson's rule gives the value $\frac{h}{3}(-h^3+0+h^3)=0$ for $\int _{-h}^hx^3dx$, whose value is also zero. This is a bonus of the rule; it actually integrates exactly a higher degree polynomial than we would expect.

Composite Simpson Rule

 

Figure: Subdivisions for Simpson's Rule

We divide the interval [a,b] into 2n subintervals of width $h=\frac{b-a}{2n}$ and apply Simpson's Rule to the n intervals $[a,a+2h],[a+2h,a+4h],\ldots ,[a+(2n-2)h,b]$. This gives

$$\int _a^bf(x)dx\approx \frac{h}{3}\{ f(a)+f(b)+4\sum _{i=1}^nf(a+(2i-1)h)+2\sum _{i=1}^{n-1}f(a+2ih)\}.$$

Sometimes we are given data for an odd number of intervals 2n+1 say. In such cases we use Simpson for the first 2n intervals in the usual way. For the last interval, we could use the Trapezium Rule, but this would lower the accuracy. We derive a special formula which uses a function value outside the interval of integration.

Let

$$\int _h^{2h}f(x)dx\approx af(0)+bf(h)+cf(2h).$$

$$\begin{array}{lrcll} f(x)=1: & h & = & a+ & b+c\\ f(x)=x: & ... ...\\ f(x)=x^2: & \frac{7}{3}h^3 & = & & bh^2+4ch^2. \end{array}$$

This leads to the formula

$$\int _h^{2h}f(x)dx\approx \frac{h}{12}\{ -f(0)+8f(h)+5f(2h)\} .$$

ExampleEvaluate $\int _0^1\cos xdx$ using the data

$$\begin{array}{ccccccc} x & 0 & 0.2 & 0.4 & 0.6 & 0.8 & 1.0\\ ... ...& 0.98007 & 0.92106 & 0.82534 & 0.69671 & 0.54030. \end{array}$$

The Trapezium Rule gives

$$\int _0^1\cos xdx\approx 0.2(\frac{1+0.54030}{2}+0.98007+\cdots +0.69671) \approx 0.83867.$$

The correct value is 0.84147.

Simpson's Rule for [0,0.8] gives

$$\frac{0.2}{3}(1+4\times 0.98007+2\times 0.92106+4\times 0.82534+0.69671) \approx 0.71736.$$

The Trapezium Rule for [0.8,1] gives

$$\frac{0.2}{2}(0.69671+0.54030)\approx 0.12370,$$

so the approximation to the whole integral is 0.84106, which is certainly better than the answer given by using the Trapezium Rule over the whole interval.

Using the special rule derived above, the integral over [0.8,1] is

$$\frac{0.2}{12}(-0.82534+8\times 0.69671+5\times 0.54030)\approx 0.12416.$$

Combining this with the Simpson estimate over [0,0.8] we obtain for the whole integral 0.84152.end.

Error Analysis
Although we can repeat calculations with smaller subintervals, it would be useful to have some idea of the errors we incur, even if we only want to compare methods.

Trapezium rule

$$\begin{eqnarray*}\mbox{local truncation error} & = & \int _0^hf(x)dx-\frac{h}{2}... ...h^3}{4}f''(0)+\cdots \} \\ & = & -\frac{h^3}{12}f''(0)+\cdots . \end{eqnarray*}$$

Simpson's Rule

$$\begin{eqnarray*}\mbox{error} & = & \int _{-h}^hf(x)dx-\frac{h}{3}\{ f(-h)+4f(0)+f(h)\} \\ & = & -\frac{h^5}{90}f^{(4)}(0)+\cdots . \end{eqnarray*}$$

This tells us that the error decreases much more rapidly for Simpson's Rule as we reduce h. To see the overall error reduction from reducing h we must look at the composite rule. For the Trapezium case

$$\begin{eqnarray*}\int _a^bf(x)dx & = & \left( \int _a^{a+h}+\int _{a+h}^{a+2h}+\... ...k=0}^n\mbox{$''$ }f(a+kh)-\frac{h^2}{12}(b-a)\bar{f''} +\cdots , \end{eqnarray*}$$

where n has been introduced in the numerator and denominator of the last term of the second line, the double prime on the last sum means halve the first and last terms and $\bar{f''}$ is the mean of the values of f''. Thus the composite error is O(h²), one degree less than the single interval error; similarly, Simpson's Rule has composite error O(h⁴). Halving the interval therefore reduces the overall error by 4 for the Trapezium Rule and by 16 for Simpson's Rule.

Clever Tricks
Let

$$I=\int _a^bf(x)dx=T(h)-kh^2 \\$$ (1)

where $k=\frac{b-a}{12}\bar{f''}$ is assumed constant, T(h) is the composite Trapezium Rule with interval h. With h halved this becomes

$$I=T(\frac{h}{2})-\frac{kh^2}{4} \\$$ (2)

Eliminating k between equations () and (), we find

$$I=\frac{4}{3}T(\frac{h}{2})-\frac{1}{3}T(h),$$

which we would expect to give a better answer.

With $a=0,b=h,\quad T(h)=\frac{h}{2}(f_0+f_h),\quad T(\frac{h}{2})=\frac{h}{4}(f_0+2f_{h/2}+f_h)$, so that

$$\frac{4}{3}T(\frac{h}{2})-\frac{1}{3}T(h)=\frac{4}{3}\frac{h}... ...\frac{1}{3}\frac{h}{2}(f_0+f_h)=\frac{h}{6}(f_0+4f_{h/2}+f_1).$$