Submitted to HEJ Manuscript no.: ANM-010201-A

Numerical integration

We have seen that some functions cannot be integrated in terms of elementary functions. This is not as bad as you may expect, because in practice it often suffices to find an approximate value of the definite integral and there exists good methods of numerical approximation. We shall now discuss the simplest and most obvious methods called quadrature techniques. We wish to direct special attention to the fundamental fact that the meaning of an approximate calculation is not precise unless it is supplemented by an estimate of the errors occurring.

The notations and conventions used here are:

$$I = \int _{a}^{b}f(x)\,dx,$$

where $a$ and $b$ are constants $(a<b)$ and $f(x)$ is given as a continuous differentiable function. Then we know that the integral exists. We have seen in Chapter 2 that we may restrict to regular subdivisions, therefore let the interval of integration be divided into $n$ equal parts of length $h=(b-a)/n$. We denote the points of subdivisions by $x_{0}=a,x_{1}=a+h,\ldots,x_{n}=b$, and the values of the function at the points of division by $f_{0},f_{1},\ldots,f_{n}$. Similarly, the values of the function at the midpoints of the intervals by ${f_{\frac {1}{2}}}$, ${f_{\frac {3}{2}}},\ldots,{f_{\frac {2\,n - 1}{2}}}$. It is fairly natural to use the Riemann sum approximation to $I$, as we did it in section 2. The procedure subdiv(a,b,n) -- we shall use it in the next examples -- returns a regular subdivision $[a=x_{0},x_{1},\ldots,x_{n}=b]$.

EXAMPLE 6.1.

$>$ f:=x->1/x: S:= subdiv(1,2,10):

$>$ integral_plot(f,S,'riemann');

 

Figure 7:

$>$ upper(f,S);

$$.7187714032$$

$>$ lower(f,S);

$$.6687714032$$

We have a lower approximation for $I$ as $0.66877$ and an upper one, which is $0.71877$. These differ by $0.05$ so there may be an error greater then 7% if either of these values is taken as an estimate for the integral. ``OK, let us simply choose a finer subdivision and we shall get a better approximation'' would be the next idea. Observe, however, that we have to calculate the minimum (or the maximum) of the function in every subinterval $[x_{i},x_{i + 1}]$, which is rather costly and error-prone. This method is easy to apply if the integrand $f(x)$ is known to be monotonic. And the other cases? The technique mentioned above is not a very good one from practical point of view. We require better techniques to satisfy two important criteria:

1. The technique should give results of high accuracy.
2. The technique should require only a small number of function evaluations.

In most cases the Riemann-sum method does not satisfy the criteria 1 and 2 but it enables us to formulate other methods realizing better these desirable properties. The most obvious method for approximating to $I$ without calculating maximum (or minimum) values on the subintervals is directly connected with the good old stepfunctions. For the value of the stepfunction on $[x_{i},x_{i + 1}]$ we choose $f_{i}$. Then for the integral we have got an approximate expression:

$$I\approx h(f_{0}+f_{1}+\ldots+f_{n - 1}).$$

This method is called rectangle rule. (Here and hereafter the symbol $\approx$ means ``is approximately equal to''.)

EXAMPLE 6.2.

$>$ f:=x->x^(1/2)*sin(x)*cos(x): n:=6: S:=subdiv(0,2,n):

$>$ integral_plot(f,S,'rectangle');

 

Figure 8:

$>$ int_num(f,S,'rectangle');

$$.3903146893$$

We obtain a closer approximation with no greater trouble if we replace the rectangle area by a trapezoid area $h (f_{i}+f_{i + 1})/2$. For the whole integral this gives the approximate expression:

$$I\approx h (f_{1}+f_{2}+\ldots+f_{n - 1})+h (f_{0}+f_{n})/2.$$

This is the trapezoid formula, since, when the areas of the trapezoids are added, each value of the function except the first and the last occurs twice.

$>$ integral_plot(f,S,'trapezoid');

 

Figure 9:

$>$ int_num(f,S,'trapezoid');

$$.3011246599$$

The approximation becomes even better if instead of choosing the trapezoid under a chord we chose the trapezoid under the tangent to the curve at the point with the abscissa $x=x_{i}+\frac {h}{2}$. The area of this trapezoid is simply $hf_{i + \frac {1}{2}}$, and the approximation for the entire integral is

$$I\approx h (f_{\frac {1}{2}}+f_{\frac {3}{2}}+\ldots+f_{\frac {2\,n - 1}{2}})$$

which is called the tangent formula. Let us see how it is work.

$>$ integral_plot(f,S,'tangent');

 

Figure 10:

$>$ int_num(f,S,'tangent');

$$.3187318652$$

The next method is called mid-point rule. In this case we assume the function values between $x_{i}$ and $x_{i + 1}$ are constant and have the value equal to the function value at $(x_{i}+x_{i + 1})/2$.

$>$ integral_plot(f,S,'mid-point');

 

Figure 11:

$>$ int_num(f,S,'mid-point');

$$.3187318652$$

Observe, that the values of the approximating integral in the case of the tangent method and the mid-point method are the same. Is it a coincidence? Explain why not.

The next method depends on estimating the subarea of the integral at two adjacent subintervals, i.e. between the abscissa $x_{i}$ and $x_{i}+2 h = x_{i + 2}$ by considering the upper boundary to be no longer a straight line but a parabola. To be more precise, the parabola which passes through the three points of the curve with abscissa ${x_{i}}$, ${x_{i + 1}}$, ${x_{i + 2}}$. The equation of the parabola is

$$y = f_{i} + \frac {(x - {x_{i}})\,({f_{i + 1}} - {f_{i}})}{h}... ...})\,(x-{x_{i}}-h)\,({f_{i+2}}-2\,{f_{i+1}}+{f_{i}})}{2\,h^{2}}.$$

Exercise 6.1. Prove the latter statement.

After a brief calculation we get the area under the parabola:

$$\frac {h\,({f_{i}} + 4\,{f_{i + 1}} + {f_{i + 2}})}{3}.$$

If we assume that $n=2m$, i.e. that $n$ is even, by the addition of the subareas we obtain the Simpson's rule:

$$I \approx \frac {h\,(4\,(\sum _{i=1}^{m}\,{f_{2\,i - 1}}) + 2\,(\sum _{i=1 }^{m - 1}\,{f_{2\,i}}) + {f_{0}} + {f_{2\,m}})}{3}.$$

$>$ integral_plot(f,S,'simpson');

 

Figure 12:

$>$ int_num(f,S,'simpson');

$$.3141154642$$

The integral of $f(x)$ can also be expressed in closed form.

$>$ Int(x^(1/2)*sin(x)*cos(x),'x'=0..2)=int(f(x),'x'=0..2);

$${\displaystyle \int _{0}^{2}} \sqrt{x}\,{\rm sin}(x)\,{\rm c... ... FresnelC}(2\, {\displaystyle \frac {\sqrt{2}}{\sqrt{\pi }}} )$$

(If you don't agree, feel free to check!) The approximate value:

$>$ evalf(int(x^(1/2)*sin(x)*cos(x),'x'=0..2));

$$.3126735508$$

It is worth to compare the results getting by the different methods:

• Rectangle rule: .3903146893
• Trapezoid rule: .3011246599
• Tangent rule: .3187318652
• Simpson's rule: .3141154642
• Maple approximation of the closed form: .3126735508

In our example the Simpson's rule seems to be the best concerning the criteria 1 and 2. Is this always true or only a coincidence again? We need an estimation of the error.

Estimations can be easily given for each of our methods, if bounds for the derivatives of the function $f(x)$ are known throughout the interval of integration. Let $M_{1},M_{2},\ldots$ be upper bounds for the absolute values of the first, second, $\ldots$ derivatives, i.e. we assume that throughout the interval $\big\vert f^{(i)}(x)\big\vert < M_{i}$ for all $i$. Then the estimation formulas are as follows:

• Rectangle rule: $\frac {{M_{1}}\,(b - a)\,h}{2}$,
• Tangent rule: $\frac {{M_{2}}\,(b - a)\,h^{2}}{24}$,
• Trapezoid rule: $\frac {{M_{2}}\,h^{3}}{12}$,
• Simpson's rule: $\frac {{M_{4}}\,h^{5}}{90}$.

From the last two estimates there also follow estimates for the entire integral. We see, that Simpson's rule has an error of much higher order in the small quantity $h$ than the others, so that where ${M_{4}}$ is not too large, it is very advantageous for practical calculations. The proofs of the estimates are quite simple using the Taylor's theorem and can be found in every textbooks on numerical methods.

There are other kinds of very important quadrature methods based on orthogonal polynomials (Gaussian quadratures). There are known (and used) methods if the integrand $f(x)$ is not finite or one of its derivatives is infinite at some point in $[a,b]$ (singular integrals). In many cases miscellaneous techniques are also available (series expansions, Laplace or Fourier transforms).

Exercise 6.2.

(a) Take a larger value for $n$, the number of subintervals, redo the numerical integrations above and check the improvements in the approximations.
(b) The numerical computations can be improved by allow a higher value to the built-in Maple variable Digits. E.g.

$>$ Digits:=20: int_num(f,S,'simpson'); Digits:=10:

$$.31411546415523086841$$

Recompute the numerical integral values of EXAMPLE 5.2 using $20$ decimal ``Digits''.
(c) Calculate

$>$ Int(exp(-x^2),'x'=0..infinity);

$${\displaystyle \int _{0}^{\infty }} e^{( - x^{2})}\,dx$$

to within $0.01$.
(d) From the formula

$>$ Pi/4 = Int(1/(1+x^2),'x'=0..1);

$${\displaystyle \frac {1}{4}} \,\pi ={\displaystyle \int _{0}^{1} } {\displaystyle \frac {1}{x^{2} + 1}} \,dx$$

calculate $\pi$

• using the trapezoid formula with $h=0.1$,
• using Simpson's rule with $h=0.1$.

(e) Calculate

$>$ Int(1/(sqrt(1+x^4)),'x'=0..1);

$${\displaystyle \int _{0}^{1}} {\displaystyle \frac {1}{\sqrt{1 + x^{4}}}} \,dx$$

numerically with an error less than $0.1$.

Exercise 6.3.

Using the trapezoid rule with $n$ subdivisions -- let us denode the integral approximation by $I_n$ -- calculate the following integrals with $n=2,4,8,16,\ldots,512$.

• $\int _0^1 e^{-x^2} dx$,
• $\int _0^1 x^{5/2} dx$,
• $\int _{-4}^4 \frac{dx}{1+x^2}$,
• $\int _0^{2\pi} \frac{dx}{2+\cos (x)}$,
• $\int _0^{\pi} e^x \cos (4x) dx$.

Exercise 6.4.

Repeat the previous exercise using Simpson's rule.

Submitted to HEJ Manuscript no.: ANM-010201-A