6.6 Numerical Integration Methods

Section 5.2 showed how to obtain exact values for definite integrals of some simple functions (low-degree polynomials) by using areas of rectangles. For functions with no closed-form antiderivative, the rectangle method typically produces an approximate value of the definite integral—the more rectangles, the better the approximation.

[r] $tikzpicture[>=latex,every node/.style=font=,scale=2.0] [fillcolor] (0,0) -- plot[domain=0:sqrt(pi),samples=200,smooth] functionsin(x*x) -- (1.7724,0) -- (0,0); [linecolor,line width=1.3pt] plot[domain=0:sqrt(pi),samples=200,smooth] functionsin(x*x); [above] at (1.253,1) y = \,(x^2); [<->,black!60,line width=1pt,anchor=base] (0,1.2) node[above] y |- (2.2,0) node[right] x node[black,shift=(0,-0.4)] at (0,0) 0 node[black,shift=(0,-0.4)] at (1.7724,0) ; [black!60,line width=1pt] (-1pt,1) -- (0pt,1); [left] at (-1pt,1) 1; tikzpicture$ For example, suppose you wanted to evaluate

\[ \int_0^{\sqrt{\pi}} \sin\,(x^2)~\dx \]

with the rectangle method. This means finding the area of the shaded region in the figure on the right. Computers have eliminated the need to do these sorts of calculations by hand. Though the rectangle method is simple to implement in a traditional programming language (e.g. via a looping construct), there are easier ways in a domain-specific language (DSL) geared toward scientific computing. One such DSL is MATLAB, or its free open-source clone Octave.^[1]

Implementing the rectangle method from scratch in Octave is a one-liner. For example, suppose you divide the interval $\ival{0}{\sqrt{\pi}}$ into 100,000 ($10^5$) subintervals of equal length, producing 100,001 (1 + 1e5 in scientific notation) equally spaced points in $\ival{0}{\sqrt{\pi}}$ (including $0$ and $\sqrt{\pi}\,$). First use the left endpoints of the $10^5$ subintervals:

octave> sum(sin(linspace(0,sqrt(pi),1+1e5)(1:end-1).^2)*sqrt(pi)/1e5)
ans = 0.8948314693913354

Now use the right endpoints:

octave> sum(sin(linspace(0,sqrt(pi),1+1e5)(2:end).^2)*sqrt(pi)/1e5)
ans = 0.8948314693913354

Finally, use the midpoints of the subintervals:

octave> sum(sin((linspace(0,sqrt(pi),1+1e5)(1:end-1)+sqrt(pi)/2e5).^2)*
        sqrt(pi)/1e5)
ans = 0.8948314695305527

The true value of the integral up to 15 decimal places is 0.894831469484145, so all three approximations are accurate to 9 decimal places. The syntax in the above commands can be explained with some examples. The following command creates 4 equally spaced points in the interval $\ival{1}{7}$ (including $x=1$ and $x=7$), thus dividing $\ival{1}{7}$ into 3 subintervals each of length $(7-1)/3 = 2$:

octave> linspace(1,7,4)
ans =

  1  3  5  7

Get all but the last number in the above list:^[2]

octave> linspace(1,7,4)(1:end-1)
ans =

  1  3  5

Now square each number in that list of numbers (the dot before the exponentiation operator ^ applies the squaring operation ^2 element-wise in the list):

octave> linspace(1,7,4)(1:end-1).^2
ans =

  1  9  25

Now take the sine of each of those squared numbers (measured in radians):

octave> sin(linspace(1,7,4)(1:end-1).^2)
ans =

  0.8414709848078965  0.4121184852417566  -0.132351750097773

Now multiply each of those numbers (the heights of the rectangles) by the width (2) of each rectangle, then add up those areas:

octave> sum(sin(linspace(1,7,4)(1:end-1).^2)*2)
ans = 2.24247543990376

Before the advent of modern computing, the rectangle method was considered inefficient, and so alternative methods were created. Two such methods are the trapezoid rule and Simpson’s rule. The idea behind both methods is to take advantage of a nonlinear function’s changing slope by using nonrectangular regions. For the trapezoid rule those regions are trapezoids, while Simpson’s rule uses quasi-rectangular regions whose top edges are parabolas, as shown in Figure 6.6.1:

tikzpicture[>=latex,every node/.style=font=] [fill=fillcolor,fill opacity=0.5,domain=0.5:2.8,samples=500,variable= ] (0.5,0) -- plot ( , 2+(exp( )*cos(4* r))/10) -- (2.8,0) -- cycle; [dashed] (0.5,0) -- (0.5,1.931); [dashed] (2.8,0) -- (2.8,2.333); [fill=blue!40,fill opacity=0.5] (1,0) -- (1,2.475) -- (1.7,2.475) -- (1.7,0); [black!60,line width=1pt,<->,anchor=base] (0,3) node[above] y |- (3.5,0) node[right] x node[black,shift=(0,-0.4)] at (0.5,0) a node[black,shift=(0,-0.4)] at (2.8,0) b node[black,shift=(0,-0.4)] at (1,0) x_i node[black,shift=(0,-0.4)] at (1.7,0) x_i+1; [linecolor,line width=1.5pt,domain=0.5:2.8,samples=500,variable= ] plot ( , 2+(exp( )*cos(4* r))/10); [black!60] (1,0) -- (1,2.475) -- (1.7,2.475) -- (1.7,0); [above] at (2.6,2.5) y=f(x); tikzpicture — Figure 6.6.1 Comparison of numerical integration methods for $\int_a^b f(x)\,\dx$

For a partition $P= \lbrace a=x_0 < x_1 < \cdots < x_{n-1} < x_n = b \rbrace$ of an interval $\ival{a}{b}$ into $n \ge 1$ subintervals of equal width $h = (b-a)/n$, let $y_i = f(x_i)$ for $i = 0$, $1$, $\ldots$, $n$. The trapezoid rule adds up the areas of trapezoids on each subinterval $\ival{x_i}{x_{i+1}}$, with the top edge being the line segment joining the points $(x_i,y_i)$ and $(x_{i+1},y_{i+1})$. The approximation formula is straightforward to derive, based on areas of trapezoids:

Theorem 6.8

Trapezoid rule:

\begin{equation*} \int_a^b f(x)~\dx ~\approx~ \frac{h}{2}\,(y_0 ~+~ 2y_1 ~+~ 2y_2 ~+~ \cdots ~+~ 2y_{n-1} ~+~ y_n) \end{equation*}

Simpson’s rule depends on pairs of neighboring subintervals: $\ival{x_0}{x_1}$ and $\ival{x_1}{x_2}$, $\ival{x_2}{x_3}$ and $\ival{x_3}{x_4}$, $\ldots$ , $\ival{x_{n-2}}{x_{n-1}}$ and $\ival{x_{n-1}}{x_n}$. Thus, $n \ge 2$ must be even. The top edge of the region over each pair $\ival{x_{i}}{x_{i+1}}$ and $\ival{x_{i+1}}{x_{i+2}}$ is the unique parabola joining the 3 points $(x_i,y_i)$, $(x_{i+1},y_{i+1})$, and $(x_{i+2},y_{i+2})$. The approximation formula is then:^[3]

Theorem 6.9

Simpson’s rule:

\begin{equation*} \int_a^b f(x)~\dx ~\approx~ \frac{h}{3}\,(y_0 ~+~ 4y_1 ~+~ 2y_2 ~+~ 4y_3 ~+~ 2y_4 ~+~ \cdots ~+~ ~+~ 2y_{n-2} ~+~ 4y_{n-1} ~+~ y_n) \end{equation*}

Example 6.32

Approximate the value of $~\displaystyle\int_0^{\sqrt{\pi}} \sin\,(x^2)~\dx~$ by using the trapezoid rule and Simpson’s rule with $n=10^5$ subintervals.

Solution: Since $x_0 = 0$ and $x_n = \sqrt{\pi}$, then $y_0 = \sin\,(x_0^2) = \sin\,0 = 0$ and $y_n = \sin\,(x_n^2) = \sin\,\pi = 0$. Thus, $y_0$ and $y_n$ contribute nothing to the summation formulas for both rules. In particular the trapezoid rule approximation becomes

\[ \int_0^{\sqrt{\pi}} \sin\,(x^2)~\dx ~\approx~ \frac{h}{2}\,(0 ~+~ 2y_1 ~+~ 2y_2 ~+~ \cdots ~+~ 2y_{n-1} ~+~ 0) ~=~ h\,\cdot\,\left(\sum_{k=1}^{n-1} y_k\right) \]

which is simple to implement in Octave:

octave> x = linspace(0,sqrt(pi),1+1e5);
octave> h = sqrt(pi)/1e5;
octave> h*sum(sin(x(2:end-1).^2))
ans = 0.8948314693913405

Likewise, the Simpson’s rule approximation becomes

\[ \int_0^{\sqrt{\pi}} \sin\,(x^2)~\dx ~\approx~ \frac{h}{3}\,(4y_1 ~+~ 2y_2 ~+~ 4y_3 ~+~ 2y_4 ~+~ \cdots ~+~ ~+~ 2y_{n-2} ~+~ 4y_{n-1}) \]

which can be implemented easily by using Octave’s powerful indexing features:

octave> (h/3)*(4*sum(sin(x(2:2:end-1).^2)) + 2*sum(sin(x(3:2:end-1).^2)))
ans = 0.8948314694841457

In the above command, the statement x(3:2:end-1) allows you to skip every other element in the list x after position 3, by moving up the list in increments of 2 positions all the way to the next-to-last position in the list (end-1). Similarly for x(2:2:end-1), which starts at position 2 and then moves up in increments of 2.

Note that Simpson’s rule gives essentially the true value in this case, and the value from the trapezoid rule is virtually the same as the value produced by the built-in trapz function in Octave/MATLAB:

octave> trapz(x,sin(x.^2))
ans = 0.8948314693913402

In general you are better off using these sorts of built-in functions instead of implementing your own.

Typically Simpson’s rule is slightly more efficient than the trapezoid rule, which is slightly more efficient than the rectangle method. However, in the above examples all the approximations were accurate to at least 9 decimal places (equivalent to getting the distance between Detroit and Chicago correct within the thickness of a toothpick). The running time of each calculation was only a few thousandths of a second. Modern computing has generally made the efficiency differences negligible. Notice that the approximations in the rectangle method, the trapezoid rule and Simpson’s rule can all be written as linear combinations of function values $f(a_i)$ multiplied by “weights” $w_i$:

\[ \int_a^b f(x)~\dx ~\approx~ \sum_{i=0}^{n} w_i\,f(a_i) \]

For example, the weights in Simpson’s rule are $w_i = \frac{h}{3}$, $\frac{2h}{3}$, or $\frac{4h}{3}$, depending on the points $a_i$ in the interval $\ival{a}{b}$. The method of Gaussian quadrature transforms an integral over any interval $\ival{a}{b}$ into an integral over the specific interval $\ival{-1}{1}$ and then uses a standard set of points in $\ival{-1}{1}$ and known weights for those points:^[4]

Theorem 6.10

Gaussian quadrature: Transform the integral $\int_a^b f(x)\,\dx$ into an integral over $\ival{-1}{1}$ by means of the substitution $u = \frac{1}{b-a}(2x - a - b)$, so that $x = \frac{b-a}{2} u + \frac{a+b}{2}$ and $\dx = \frac{b-a}{2}\du$. Then

\[ \int_a^b f(x)~\dx ~=~ \frac{b-a}{2}\,\int_{-1}^1 g(u)~\du ~\approx~ \frac{b-a}{2}\,\sum_{i=1}^n w_i\,g(a_i) \]

where $g(u) = f\left(\frac{b-a}{2}\,u + \frac{a+b}{2}\right)$, with the points $a_1$, $\ldots$, $a_n$ and weights $w_1$, $\ldots$, $w_n$ given in Table 6.6 below for any choice of $2 \le n \le 10$ points in $\ival{-1}{1}$.

Table 6.1 Table of Gaussian quadrature points and weights

$n$	$a_1$, $\ldots$, $a_n$	$w_1$, $\ldots$, $w_n$
$2$	$\pm 0.577350$	$1$
$3$	$\hphantom{\pm} 0$	$8/9$
	$\pm 0.774597$	$5/9$
$4$	$\pm 0.339981$	$0.652145$
	$\pm 0.861136$	$0.347855$
$5$	$\hphantom{\pm} 0$	$0.568889$
	$\pm 0.538469$	$0.478629$
	$\pm 0.906180$	$0.236927$
$6$	$\pm 0.238619$	$0.467914$
	$\pm 0.661209$	$0.360762$
	$\pm 0.932470$	$0.171324$
$7$	$\hphantom{\pm} 0$	$0.417959$
	$\pm 0.405845$	$0.381830$
	$\pm 0.741531$	$0.279705$
	$\pm 0.949108$	$0.129485$

$n$	$a_1$, $\ldots$, $a_n$	$w_1$, $\ldots$, $w_n$
$8$	$\pm 0.183435$	$0.362684$
	$\pm 0.525532$	$0.313707$
	$\pm 0.796666$	$0.222381$
	$\pm 0.960290$	$0.101229$
$9$	$\hphantom{\pm} 0$	$0.330239$
	$\pm 0.324253$	$0.312347$
	$\pm 0.613371$	$0.260611$
	$\pm 0.836031$	$0.180648$
	$\pm 0.968160$	$0.081274$
$10$	$\pm 0.148874$	$0.295524$
	$\pm 0.433395$	$0.269267$
	$\pm 0.679410$	$0.219086$
	$\pm 0.865063$	$0.149451$
	$\pm 0.973907$	$0.066671$

Example 6.33

Approximate the value of $~\displaystyle\int_0^2 \dfrac{\dx}{1 + x^3}~$ by using Gaussian quadrature with $n=4$ points.

Solution: For $a=0$ and $b=2$, use the substitution $u = \frac{1}{b-a}(2x - a - b) = x-1$, so that $x=u+1$ and $\dx = \du$. Thus, $g(u) = f(u+1) = \frac{1}{1 + (u+1)^3}$. Using $n=4$ in Table 6.6, the points $a_i$ and weights $w_i$ are

\begin{alignat*}{3} a_1 ~&=~ -0.339981 \quad\quad & w_1 ~&=~ 0.652145\\ a_2 ~&=~ 0.339981 \quad\quad & w_2 ~&=~ 0.652145\\ a_3 ~&=~ -0.861136 \quad\quad & w_3 ~&=~ 0.347855\\ a_4 ~&=~ 0.861136 \quad\quad & w_4 ~&=~ 0.347855 \end{alignat*}

and so

\begin{align*} \int_0^2 \frac{\dx}{1 + x^3} ~&=~ \frac{2-0}{2}\,\int_{-1}^1 g(u)~\du ~=~ \int_{-1}^1 \frac{\du}{1 + (u+1)^3} ~\approx~ \sum_{i=1}^4 w_i\,g(a_i) ~=~ \sum_{i=1}^4 w_i\,\cdot\,\frac{1}{1 + (a_i+1)^3}\\[6pt] &\approx~ \frac{0.652145}{1 + (-0.339981+1)^3} ~+~ \frac{0.652145}{1 + (0.339981+1)^3} ~+~ \frac{0.347855}{1 + (-0.861136+1)^3} ~+~ \frac{0.347855}{1 + (0.861136+1)^3}\\[4pt] &\approx~ 1.091621 \end{align*}

The true value of the integral to six decimal places is 1.090002.
Using more points (e.g. $n=7$) is easy to implement in Octave, using element-wise operations on arrays:

octave> a = [ 0 -0.405845 0.405845 -0.741531 0.741531 -0.949108 0.949108 ];
octave> w = [ 0.417959 0.381830 0.381830 0.279705 0.279705 0.129485 0.129485 ];
octave> sum(w./(1 + (a+1).^3))
ans = 1.090016688064804

Gaussian quadrature can be applied to improper integrals. For example,

\[ \int_0^{\infty} f(x)\,e^{-x}\,\dx ~\approx~ \sum_{i=1}^n w_i\,f(a_i) \]

using the points $a_i$ and weights $w_i$ in Table 6.6^[5] for $n=3, 4$, or 5 points in $\lival{0}{\infty}$:

Table 6.2 Table of Gaussian quadrature points and weights for $\int_0^{\infty} f(x)\,e^{-x}\,\dx$

$n$	$a_1$, $\ldots$, $a_n$	$w_1$, $\ldots$, $w_n$
$3$	$0.415775$	$0.711093$
	$2.294280$	$0.278518$
	$6.289945$	$0.010389$

$n$	$a_1$, $\ldots$, $a_n$	$w_1$, $\ldots$, $w_n$
$4$	$0.322548$	$0.603154$
	$1.745761$	$0.357419$
	$4.536620$	$0.038888$
	$9.395071$	$0.000539$

$n$	$a_1$, $\ldots$, $a_n$	$w_1$, $\ldots$, $w_n$
$5$	$0.263560$	$0.521756$
	$1.413403$	$0.398667$
	$3.596426$	$0.075942$
	$7.085810$	$0.003612$
	$12.640801$	$0.000023$

Example 6.34

Approximate the value of $~\displaystyle\int_0^{\infty} x^5\,e^{-x}\,\dx~$ by using Gaussian quadrature with $n=3$ points in Table 6.6.

Solution: For $n=3$, Table 6.6 gives $a_1=0.415775$, $a_2=2.294280$, $a_3=6.289945$, and $w_1=0.711093$, $w_2=0.278518$, $w_3=0.010389$. Then for $f(x)=x^5$,

\begin{align*} \int_0^{\infty} x^5\,e^{-x}\,\dx ~&\approx~ \sum_{i=1}^n w_i\,f(a_i) ~=~ \sum_{i=1}^n w_i\,a_i^5\\[4pt] &\approx~ 0.711093\,(0.415775)^5 ~+~ 0.278518\,(2.294280)^5 ~+~ 0.010389\,(6.289945)^5\\ &\approx~ 119.9974709727211 \end{align*}

The true value is $\Gamma\,(6) = 5! = 120$.
Note: The points $a_i$ in Table 6.6 are the roots of the Laguerre polynomials of degree $n$.

Exercises

A simple pendulum of length $l$ swings through an angle of $90\Degrees$ on either side of the vertical with period $P$, given by

\[ P ~=~ 4\,\sqrt{\dfrac{l}{g}}\,\int_0^{\pi/2} \frac{\dtheta}{\sqrt{1 \;-\; 0.5\,\sin^2 \theta}} \]

where $g = 9.8$ m/s² is the acceleration due to gravity. Use the rectangle method (with left endpoints), the trapezoid rule, and Simpson’s rule to write $P$ as a constant multiple of $\sqrt{l/g}$. Preferably, use a computer and $n= 10^5$ subintervals of equal width (or $n=10$ subintervals if calculating by hand).
Repeat Exercise 1 using Gaussian quadrature with $n=5$ points.
Approximate the value of $~\displaystyle\int_0^{\sqrt{\pi}} \sin\,(x^2)~\dx~$ by using Gaussian quadrature with $n=7$ points. 2
Repeat Exercise 3 with $n=9$ points.
Repeat Exercise 3 with $n=10$ points.
The points $a_i$ in Table 6.6 for Gaussian quadrature are the roots of the Legendre polynomials $P_n(x)$, with $P_0(x) = 1$, $P_1(x) = x$, and $P_n(x)$ defined for integers $n \ge 2$ by the recursion formula

\[ n\,P_{n}(x) ~=~ (2n-1)\,x\,P_{n-1}(x) ~-~ (n-1)\,P_{n-2}\,(x) \]
1. Write out $P_n(x)$ explicitly in standard polynomial form for $n = 2, 3, 4, 5$.
2. Verify that the roots of $P_n(x)$ match the $n$ points $a_1$, $\ldots$, $a_n$ in Table 6.6 for $n = 2, 3, 4, 5$.
3. With no calculations, explain why $~\int_{-1}^1 P_2(x)\,P_3(x)\,\dx = \int_{-1}^1 P_3(x)\,P_4(x)\,\dx = \int_{-1}^1 P_4(x)\,P_5(x)\,\dx = 0$.
4. For $n=0, 1, 2, 3$, verify that
  
  \[ \int_{-1}^1 P_n^2(x)~\dx ~=~ \frac{2}{2n+1} ~. \]
Repeat Example 6.34 with $n=4$ points.
Repeat Example 6.34 with $n=5$ points.
Use Table 6.6 to approximate the value of $~\displaystyle\int_0^{\infty} \ln\,(1+x)\,e^{-x}\,\dx~$ with $n=5$ points.

Octave is freely available at https://www.gnu.org/software/octave/ ↩

MATLAB would require you to use two steps:

MATLAB>> x = linspace(1,7,4);
MATLAB>> x(1:end-1)

↩

For a full derivation of both formulas, see pp.144-149 in Hornbeck, R.W., Numerical Methods, New York: Quantum Publishers, Inc., 1975. ↩
The details are beyond the scope of this book. See Chapter 4 in Ralston, A. and P. Rabinowitz, A First Course in Numerical Analysis, 2nd ed., New York: McGraw-Hill, Inc., 1978. See also Table 1 in Stroud, A.H. and D. Secrest, Gaussian Quadrature Formulas, Englewood Cliffs, NJ: Prentice-Hall, Inc., 1966. ↩
See Table 6 in Stroud, A.H. and D. Secrest, Gaussian Quadrature Formulas. ↩

\(n\)	\(a_1\), \(\ldots\), \(a_n\)	\(w_1\), \(\ldots\), \(w_n\)
\(2\)	\(\pm 0.577350\)	\(1\)
\(3\)	\(\hphantom{\pm} 0\)	\(8/9\)
	\(\pm 0.774597\)	\(5/9\)
\(4\)	\(\pm 0.339981\)	\(0.652145\)
	\(\pm 0.861136\)	\(0.347855\)
\(5\)	\(\hphantom{\pm} 0\)	\(0.568889\)
	\(\pm 0.538469\)	\(0.478629\)
	\(\pm 0.906180\)	\(0.236927\)
\(6\)	\(\pm 0.238619\)	\(0.467914\)
	\(\pm 0.661209\)	\(0.360762\)
	\(\pm 0.932470\)	\(0.171324\)
\(7\)	\(\hphantom{\pm} 0\)	\(0.417959\)
	\(\pm 0.405845\)	\(0.381830\)
	\(\pm 0.741531\)	\(0.279705\)
	\(\pm 0.949108\)	\(0.129485\)

\(n\)	\(a_1\), \(\ldots\), \(a_n\)	\(w_1\), \(\ldots\), \(w_n\)
\(8\)	\(\pm 0.183435\)	\(0.362684\)
	\(\pm 0.525532\)	\(0.313707\)
	\(\pm 0.796666\)	\(0.222381\)
	\(\pm 0.960290\)	\(0.101229\)
\(9\)	\(\hphantom{\pm} 0\)	\(0.330239\)
	\(\pm 0.324253\)	\(0.312347\)
	\(\pm 0.613371\)	\(0.260611\)
	\(\pm 0.836031\)	\(0.180648\)
	\(\pm 0.968160\)	\(0.081274\)
\(10\)	\(\pm 0.148874\)	\(0.295524\)
	\(\pm 0.433395\)	\(0.269267\)
	\(\pm 0.679410\)	\(0.219086\)
	\(\pm 0.865063\)	\(0.149451\)
	\(\pm 0.973907\)	\(0.066671\)

\(n\)	\(a_1\), \(\ldots\), \(a_n\)	\(w_1\), \(\ldots\), \(w_n\)
\(3\)	\(0.415775\)	\(0.711093\)
	\(2.294280\)	\(0.278518\)
	\(6.289945\)	\(0.010389\)

\(n\)	\(a_1\), \(\ldots\), \(a_n\)	\(w_1\), \(\ldots\), \(w_n\)
\(4\)	\(0.322548\)	\(0.603154\)
	\(1.745761\)	\(0.357419\)
	\(4.536620\)	\(0.038888\)
	\(9.395071\)	\(0.000539\)

\(n\)	\(a_1\), \(\ldots\), \(a_n\)	\(w_1\), \(\ldots\), \(w_n\)
\(5\)	\(0.263560\)	\(0.521756\)
	\(1.413403\)	\(0.398667\)
	\(3.596426\)	\(0.075942\)
	\(7.085810\)	\(0.003612\)
	\(12.640801\)	\(0.000023\)