Prev Next

Michael Oberguggenberger and Alexander OstermannAnalysis for Computer ScientistsUndergraduate Topics in Computer Sciencehttps://doi.org/10.1007/978-3-319-91155-7_13

13. Numerical Integration

Michael Oberguggenberger¹ and Alexander Ostermann¹

(1)

University of Innsbruck, Innsbruck, Austria

Michael Oberguggenberger (Corresponding author)

Email: michael.oberguggenberger@uibk.ac.at

Alexander Ostermann

Email: alexander.ostermann@uibk.ac.at

The fundamental theorem of calculus suggests the following approach to the calculation of definite integrals: one determines an antiderivative F of the integrand f and computes from that the value of the integral

$\int ^b_a \! f(x) \, \mathrm{d}x = F(b) - F(a).$

In practice, however, it is difficult and often even impossible to find an antiderivative F as a combination of elementary functions. Apart from that, antiderivatives can also be fairly complex, as the example $\int x^{100} \sin x \,\mathrm{d}x$ shows. Finally, in concrete applications the integrand is often given numerically and not by an explicit formula. In all these cases one reverts to numerical methods. In this chapter the basic concepts of numerical integration (quadrature formulas and their order) are introduced and explained. By means of instructive examples we analyse the achievable accuracy for the Gaussian quadrature formulas and the required computational effort.

13.1 Quadrature Formulas

For the numerical computation of $\int ^b_a \!f(x) \, \mathrm{d}x$ we first split the interval of integration [a, b] into subintervals with grid points $a \,{=}\, x_0< x_1< x_2< \ldots< x_{N-1} < x_N \,{=}\, b$ , see Fig. 13.1. From the additivity of the integral (Proposition 11.10 (d)) we get

$\int ^b_a \! f(x) \,\mathrm{d}x = \sum ^{N-1}_{j=0} \int ^{x_{j+1}}_{x_j} \!\!f(x) \,\mathrm{d}x.$

Hence it is sufficient to find an approximation formula for a (small) subinterval of length $h_j = x_{j+1}-x_j$ . One example of such a formula is the trapezoidal rule through which the area under the graph of a function is approximated by the area of the corresponding trapezoid (Fig. 13.2)

$\int _{x_j}^{x_{j+1}} \!\!f(x) \, \mathrm{d}x \ \approx \ h_j\, \frac{1}{2} \Bigl (f(x_j) + f(x_{j+1})\Bigr ).$

For the derivation and analysis of such approximation formulas it is useful to carry out a transformation onto the interval [0, 1]. By setting $x = x_j + \tau h_j$ one obtains from $\mathrm{d}x = h_j \, \mathrm{d}\tau$ that

$\int ^{x_{j+1}}_{x_j}\! f(x)\, \mathrm{d}x = \int ^1_0 \!f(x_j + \tau h_j) h_j \, \mathrm{d}\tau = h_j \int ^1_0 \!g(\tau ) \, \mathrm{d}\tau$

images/215236_2_En_13_Chapter/215236_2_En_13_Fig1_HTML.gif — Fig. 13.1
Partition of the interval of integration into subintervals

images/215236_2_En_13_Chapter/215236_2_En_13_Fig2_HTML.gif — Fig. 13.2
Trapezoidal rule

with $g(\tau ) = f(x_j + \tau h_j)$ . Thus it is sufficient to find approximation formulas for $\int ^1_0 g(\tau )\, \mathrm{d}\tau$ . The trapezoidal rule in this case is

$\int ^1_0 \!g(\tau ) \, \mathrm{d}\tau \ \approx \ \frac{1}{2} \Bigl (g(0) + g(1)\Bigr ).$

Obviously, it is exact if $g(\tau )$ is a polynomial of degree 0 or 1.

In order to obtain a more accurate formula, we demand that quadratic polynomials are integrated exactly as well. For the moment let

$g(\tau ) = \alpha + \beta \tau + \gamma \tau ^2$

be a general polynomial of degree 2. Due to $g(0)=\alpha$ , $g\big (\frac{1}{2}\big )=\alpha +\frac{1}{2}\beta +\frac{1}{4}\gamma$ and $g(1)=\alpha +\beta +\gamma$ we get by a short calculation

$\int ^1_0 \bigl (\alpha + \beta \tau + \gamma \tau ^2\bigr )\, \mathrm{d}\tau \ = \ \alpha + \frac{1}{2} \beta + \frac{1}{3} \gamma \ = \ \frac{1}{6} \Big (g(0) + 4g\bigl (\tfrac{1}{2}\bigr ) + g(1)\Big ).$

The corresponding approximation formula for general g reads

$\int ^1_0 g(\tau ) \, \mathrm{d}\tau \ \approx \ \frac{1}{6} \Big (g(0) + 4g\big (\tfrac{1}{2}\big ) + g(1)\Big ).$

By construction, it is exact for polynomials of degree less than or equal to 2; it is called Simpson’s rule.¹

The special forms of the trapezoidal and of Simpson’s rule motivate the following definition.

Definition 13.1

The approximation formula

$\int ^1_0 \! g(\tau ) \, \mathrm{d}\tau \ \approx \ \sum ^s_{i=1} b_i \, g(c_i)$

is called a quadrature formula. The numbers $b_1,\ldots , b_s$ are called weights, and the numbers $c_1,\ldots , c_s$ are called nodes of the quadrature formula; the integer s is called the number of stages.

A quadrature formula is determined by the specification of the weights and nodes. Thus we denote a quadrature formula by $\{(b_i, c_i),\ i=1,\ldots , s\}$ for short. Without loss of generality the weights $$b_i$$ are not zero, and the nodes are pairwise different $(c_i \ne c_k$ for $i \ne k)$ .

Example 13.2

(a) The trapezoidal rule has $$s = 2$$

stages and is given by

$b_1 = b_2 = \frac{1}{2}, \quad c_1 = 0, \quad c_2 = 1.$

(b) Simpson’s rule has $$s = 3$$

stages and is given by

$b_1 = \frac{1}{6}, \quad b_2 = \frac{2}{3}, \quad b_3 = \frac{1}{6}, \quad c_1 = 0, \quad c_2 = \frac{1}{2}, \quad c_3 = 1.$

In order to compute the original integral $\int ^b_a f(x) \, \mathrm{d}x$ by quadrature formulas, one has to reverse the transformation from f to g. Due to $g (\tau ) = f(x_j + \tau h_j)$ one obtains

$\int ^{x_{j+1}}_{x_j} \!f(x)\, \mathrm{d}x = h_j \int ^1_0 \!g(\tau ) \, \mathrm{d}t \ \approx \ h_j \sum ^s_{i=1} b_i g(c_i) = h_j \sum ^s_{i=1} b_i f(x_j + c_i h_j),$

and thus the approximation formula

$\int ^b_a \!f(x) \, \mathrm{d}x \ = \ \sum ^{N-1}_{j=0} \int ^{x_{j+1}}_{x_j} \!f(x) \, \mathrm{d}x \ \approx \ \sum ^{N-1}_{j=0} h_j \sum ^s_{i=1} b_i f(x_j + c_i h_j).$

We now look for quadrature formulas that are as accurate as possible. Since the integrand is typically well approximated by Taylor polynomials on small intervals, a good quadrature formula is characterised by the property that it integrates exactly as many polynomials as possible. This idea motivates the following definition.

Definition 13.3

(Order) The quadrature formula $\{(b_i, c_i),\ i = 1,\ldots , s\}$ has order p if all polynomials g of degree less or equal to $$p-1$$

are integrated exactly by the quadrature formula; i.e.,

$\int ^1_0 g(\tau )\, \mathrm{d}\tau = \sum ^s_{i=1} b_i\, g(c_i)$

for all polynomials g of degree smaller than or equal to $$p-1$$

Example 13.4

(a) The trapezoidal rule has order 2.

(b) Simpson’s rule has (by construction) at least order 3.

The following proposition yields an algebraic characterisation of the order of quadrature formulas.

Proposition 13.5

A quadrature formula $\{(b_i, c_i),\ i=1,\ldots , s\}$ has order p if and only if

$\sum ^s_{i=1} b_i \, c_i^{q-1} = \frac{1}{q} \qquad \text {for} \quad 1 \le q \le p.$

Proof

One uses the fact that a polynomial g of degree $$p-1$$

$g(\tau ) = \alpha _0 + \alpha _1 \tau + \ldots +\alpha _{p-1} \, \tau ^{p-1}$

is a linear combination of monomials, and that both integration and application of a quadrature formula are linear processes. Thus it is sufficient to prove the result for the monomials

$g(\tau ) = \tau ^{q-1}, \qquad 1 \le q \le p.$

The proposition now follows directly from the identity

$\frac{1}{q} = \int ^1_0 \!\tau ^{q-1} \, \mathrm{d}\tau = \sum ^s_{i=1} b_i \, g(c_i) = \sum ^s_{i=1} b_i c_i^{q-1}.{\square }$

The conditions of the proposition

$\begin{aligned} b_1 + b_2 + \ldots +b_s&= 1\\[2pt] b_1 c_1 + b_2 c_2 + \ldots + b_s c_s&= \tfrac{1}{2}\\[2pt] b_1 c_1^2 + b_2 c_2^2 + \ldots + b_s c_s^2&= \tfrac{1}{3}\\ \vdots&\\ b_1 c_1^{p-1} + b_2 c_2^{p-1} + \ldots + b_s c_s^{p-1}&=\tfrac{1}{p} \end{aligned}$

are called order conditions of order p. If s nodes $c_1,\ldots , c_s$ are given then the order conditions form a linear system of equations for the unknown weights $$b_i$$

. If the nodes are pairwise different then the weights can be determined uniquely from that. This shows that for s different nodes there always exists a unique quadrature formula of order $p\ge s$ .

Example 13.6

We determine once more the order of Simpson’s rule. Due to

$\begin{aligned} b_1 + b_2 + b_3&= \tfrac{1}{6}+\tfrac{2}{3}+\tfrac{1}{6} =1\\[2pt] b_1 c_1 + b_2 c_2 + b_3 c_3&= \tfrac{2}{3}\cdot \tfrac{1}{2}+\tfrac{1}{6} =\tfrac{1}{2}\\[2pt] b_1 c_1^2 + b_2 c_2^2 + b_3 c_3^2&= \tfrac{2}{3}\cdot \tfrac{1}{4}+\tfrac{1}{6} = \tfrac{1}{3} \end{aligned}$

its order is at least 3 (as we already know from the construction). However, additionally

$b_1 c_1^3 + b_2 c_2^3 + b_3 c_3^3 = \tfrac{4}{6} \cdot \tfrac{1}{8} + \tfrac{1}{6} = \tfrac{3}{12} = \tfrac{1}{4},$

i.e., Simpson’s rule even has order 4.

The best quadrature formulas (high accuracy with little computational effort) are the Gaussian quadrature formulas. For that we state the following result whose proof can be found in [23, Chap. 10, Corollary 10.1].

Proposition 13.7

There is no quadrature formula with s stages of order $$p > 2s$$ . On the other hand, for every $s \in \mathbb N$ there exists a (unique) quadrature formula of order $$p = 2s$$ . This formula is called s-stage Gaussian quadrature formula.

The Gaussian quadrature formulas for $s \le 3$ are

$\begin{aligned} s = 1:\quad&c_1 = \frac{1}{2}, \quad b_1 = 1, \quad \text {order 2} \ \ \text {(midpoint rule);}\\[1mm] s = 2:\quad&c_1 = \frac{1}{2}-\frac{\sqrt{3}}{6}, \quad c_2 = \frac{1}{2} + \frac{\sqrt{3}}{6},\quad \ b_1 = b_2 = \frac{1}{2}, \quad \text {order 4;}\\[1mm] s = 3:\quad&c_1 = \frac{1}{2} - \frac{\sqrt{15}}{10}, \quad c_2 = \frac{1}{2}, \quad c_3 = \frac{1}{2} + \frac{\sqrt{15}}{10},\\&b_1 = \frac{5}{18}, \quad b_2 = \frac{8}{18}, \quad b_3 = \frac{5}{18}, \quad \text {order 6}. \end{aligned}$

13.2 Accuracy and Efficiency

In the following numerical experiment the accuracy of quadrature formulas will be illustrated. With the help of the Gaussian quadrature formulas of order 2, 4 and 6 we compute the two integrals

$\int _0^3 \!\cos x \,\mathrm{d}x =\sin 3\qquad \text {and}\qquad \int _0^1\! x^{5/2}\,\mathrm{d}x = \frac{2}{7}.$

For that we choose equidistant grid points

$x_j = a+jh,\qquad j=0,\ldots , N$

with

and $N=1,2,4,8,16,\ldots , 512$ . Finally, we plot the costs of the calculation as a function of the achieved accuracy in a double-logarithmic diagram.

A measure for the computational cost of a quadrature formula is the number of required function evaluations, abbreviated by fe. For an s-stage quadrature formula, it is the number

$\mathtt{fe} = s\cdot N.$

The achieved accuracy err is the absolute value of the error. The according results are presented in Fig. 13.3. One makes the following observations:

images/215236_2_En_13_Chapter/215236_2_En_13_Fig3_HTML.gif — Fig. 13.3
Accuracy-cost-diagram of the Gaussian quadrature formulas. The crosses are the results of the one-stage Gaussian method of order 2, the squares the ones of the two-stage method of order 4 and the circles the ones of the three-stage method of order 6

(a)

The curves are straight lines (as long as one does not get into the range of rounding errors, like with the three-stage method in the left picture).
(b)

In the left picture the straight lines have slope , where p is the order of the quadrature formula. In the right picture this is only true for the method of order 2, and the other two methods result in straight lines with slope .
(c)

For given costs the formulas of higher order are more accurate.

In order to understand this behaviour, we expand the integrand into a Taylor series. On the subinterval $[\alpha , \alpha +h]$ of length h we obtain

$f(\alpha + \tau h) = \sum _{q=0}^{p-1} \frac{h^q}{q!} f^{(q)}(\alpha ) \tau ^q + {\mathcal O}(h^p).$

Since a quadrature formula of order p integrates polynomials of degree less than or equal to $$p-1$$

exactly, the Taylor polynomial of f of degree $$p-1$$

is being integrated exactly. The error of the quadrature formula on this subinterval is proportional to the length of the interval times the size of the remainder term of the integrand, so

$h\cdot {\mathcal O}(h^p) = {\mathcal O}(h^{p+1}).$

In total we have N subintervals; hence, the total error of the quadrature formula is

$N\cdot {\mathcal O}(h^{p+1}) = Nh \cdot {\mathcal O}(h^p) = (b-a)\cdot {\mathcal O}(h^p) = {\mathcal O}(h^p).$

Thus we have shown that (for small h) the error err behaves like

$\mathtt{err} \approx c_1\cdot h^p.$

Since furthermore

$\mathtt{fe} = sN = s\cdot Nh\cdot h^{-1} = s\cdot (b-a)\cdot h^{-1} = c_2\cdot h^{-1}$

holds true, we obtain

$\log (\mathtt{fe}) = \log c_2 -\log h\qquad \text {and}\qquad \log (\mathtt{err}) \approx \log c_1 + p\cdot \log h,$

so altogether

$\log (\mathtt{fe}) \approx c_3 - \frac{1}{p}\cdot \log (\mathtt{err}).$

This explains why straight lines with slope $$-1/p$$

appear in the left picture.

In the right picture it has to be noted that the second derivative of the integrand is discontinuous at 0. Hence the above considerations with the Taylor series are not valid anymore. The quadrature formula also detects this discontinuity of the high derivatives and reacts with a so-called order reduction ; i.e., the methods show a lower order (in our case $$p=7/2$$ ).

Experiment 13.8

Compute the integrals

$\int _0^3 \!\sqrt{x}\,\mathrm{d}x \quad \text {and}\quad \int _1^2\frac{\mathrm{d}x}{x}\,$

using the Gaussian quadrature formulas and generate an accuracy-cost-diagram. For that purpose modify the programs mat13_1.m, mat13_2.m, mat13_3.m, mat13_4.m and mat13_5.m with which Fig. 13.3 was produced.

Commercial programs for numerical integration determine the grid points adaptively based on automatic error estimates. The user can usually specify the desired accuracy. In MATLAB the routines quad.m and quadl.m serve this purpose.

13.3 Exercises

1.

For the calculation of $\int _0^1 x^{100} \sin x \,\mathrm{d}x$ first determine an antiderivative F of the integrand f using maple. Then evaluate $$F(1)-F(0)$$ to 10, 50, 100, 200 and 400 digits and explain the surprising results.

2.

Determine the order of the quadrature formula given by

$b_1 = b_4 = \frac{1}{8}, \quad b_2 = b_3 = \frac{3}{8},\quad c_1 = 0, \quad c_2 = \frac{1}{3}, \quad c_3 = \frac{2}{3},\quad c_4 = 1.$

3.

Determine the unique quadrature formula of order 3 with the nodes

$c_1 = \frac{1}{3}, \quad c_2 = \frac{2}{3}, \quad c_3 = 1.$

4.

Determine the unique quadrature formula with the nodes

$c_1 = \frac{1}{4}, \quad c_2 = \frac{1}{2}, \quad c_3 = \frac{3}{4}.$

Which order does it have?

5.

Familiarise yourself with the MATLAB programs quad.m and quadl.m for the computation of definite integrals and test the programs for

$\int ^1_0 \!\mathrm{e}^{-x^2}\,\mathrm{d}x \qquad \text {and} \qquad \int ^1_0 \!\root 3 \of {x} \, \mathrm{d}x.$

6.

Justify the formulas

$\pi = 4 \int ^1_0 \frac{\mathrm{d}x}{1 + x^2} \qquad \text {and} \qquad \pi = 4 \int ^1_0 \!\sqrt{1- x^2} \, \mathrm{d}x$

and use them to calculate $\pi$ by numerical integration. To do so divide the interval [0, 1] into N equally large parts $(N = 10,100,\ldots )$ and use Simpson’s rule on those subintervals. Why are the results obtained with the first formula always more accurate?

7.

Write a MATLAB program that allows you to evaluate the integral of any given (continuous) function on a given interval [a, b], both by the trapezoidal rule and by Simpson’s rule. Use your program to numerically answering the questions of Exercises 7–9 from Sect. 11.4 and Exercise 5 from Sect. 12.4.

8.

Use your program from Exercise 7 to produce tables (for $$x = 0$$

in steps of 0.5) of some higher transcendental functions:

(a): the Gaussian error function

$\text {Erf} (x) = \frac{2}{\sqrt{\pi }}\int _0^x \mathrm{e}^{-y^2} \,{\mathrm{d}}y,$
(b): the sine integral

$\mathcal {S}i (x) = \int _0^x \frac{\sin y}{y}\, \,{\mathrm{d}}y,$
(c): the Fresnel integral

$\mathcal {S}(x) = \int _0^x \sin \left( \displaystyle \frac{\pi }{2}\, y^2\right) \,{\mathrm{d}}y.$

9.

(Experimental determination of expectation values) The family of standard beta distributions on the interval [0, 1] is defined through the probability densities

$f(x;r, s) = \frac{1}{B(r,s)}\, x^{r-1}(1-x)^{s-1}, \quad 0 \le x \le 1,$

where

. Here $B(r, s) = \int _0^1 y^{r-1}(1-y)^{s-1}\,\,{\mathrm{d}}y$ is the beta function, which is a higher transcendental function for non-integer values of r, s. For integer values of $r, s \ge 1$ it is given by

$B(r, s) = \frac{(r-1)!(s-1)!}{(r+s-1)!}.$

With the help of the MATLAB program quad.m, compute the expectation values $\mu (r, s) = \int _0^1 xf(x;r, s)\,{\mathrm{d}}x$ for various integer values of r and s and guess a general formula for $\mu (r, s)$ from your experimental results.

Prev Next