Skip to main content

Section 6 Numerical Integration

Subsection 6.1 Integration review

Before we begin looking at numerical calculus, it is useful to recall some of the basic notions. In particular, we'll be reexamining the introduction to calculus most students see in the second semester of the course.

Here is a question - what is \(\int_0^1 x^2 \, dx\text{?}\) Most people that have some experience with calculus will perform the following operation:

\begin{equation*} \int_0^1 x^2 \, dx = \frac{1}{3} x^3 |_0^1 = \frac{1}{3}. \end{equation*}

However, this computation both misses the point of what the integral represents and uses a technique that will largely be unavailable in practice.

As to the first point, the definite integral represents a measurement of the signed area between a graph and the \(x\)-axis. We say signed area because area above the axis and area below the axis are considered to have opposite signs. Again, a definite integral is an area. Below, we plot the region corresponding to \(\int_0^1 x^2 \, dx\text{.}\)

Now, the area of this region is certainly \(\frac{1}{3}\text{.}\) To arrive at that conclusion, we used one of the most important theorems of continuous mathematics, the fundamental theorem of calculus, which gives a connection between definite integrals, that is the signed area under a curve, with indefinite integrals, that is antiderivatives.

One salient assumption present in the fundamental theorem of calculus is the existence of an antiderivative. Unfortunately, there are many functions that do not possess a (closed form) antiderivative, including some of the most useful functions in practice. For example, the normal distribution from statistics is essentially defined by the function \(f(x) = e^{-x^2}\text{.}\) A typical problem might wish to compute an integral like \(\int_{-1}^{2} e^{-x^2} \, dx\text{,}\) which is a simple area contained under a very nice curve, as shown below.

However, the fundamental theorem cannot be used to compute the area indicated by the definite integral because the function \(e^{-x^2}\) has no closed form antiderivative. So we need to approach the area finding problem with techniques related to the definition of definite integrals, which consist of breaking the area under functions up into approximating rectangles and then making the rectangles uniformly smaller in width.

The Riemann sum that approximates the signed area under \(f\) on \([a,b]\) is given by the following formula. Let \(n\) be the number of approximating rectangles, and the width of each rectangle be \(\Delta x = \frac{b - a}{n}.\) Then we can define a partition of \([a,b]\) by \(x_0 = a\text{,}\) \(x_i = x_0 + i \Delta x\text{,}\) and \(x_n = b\text{.}\) On each subinterval \([x_i, x_{i+1}]\text{,}\) we choose a point \(x_i^*\text{.}\) Then the area under \(f\) can be approximated by the expression

\begin{equation*} \int_a^b f(x) \, dx \approx \sum_{i = 0}^{n-1} f(x_i^*) \Delta x. \end{equation*}

Those functions for which \(\lim_{n \to \infty} \sum_{i = 0}^{n-1} f(x_i^*) \Delta x\) converges are called Riemann integrable.

Note, it need not be the case that the rectangles have equal width, which is an assumption made here to simplify the presentation.

Subsection 6.2 The trapezoid rule

An immediate observation of a Riemann sum approximation for a definite integral might lead you to conclude that other shapes might provide more accurate approximations than rectangles. An easy shape to work with in this context is the trapezoid - it has a simple formula for area and allows us to avoid having to choose random points inside the subintervals. Compare the following pictures.

While the example might seem to be artificially chosen to make the trapezoids significantly more accurate than the rectangles, in fact, the vast majority of graphs of interest will look like the pictures above for small enough subintervals. This motivates the development of the trapezoid rule for approximating a definite integral.

Recall that the area of a trapezoid with height \(h\) and base widths \(b_1, b_2\) is given by the formula

\begin{equation*} A = \frac{b_1 + b_2}{2} h. \end{equation*}

Suppose that \(f\) is a function defined on the interval \(I = [a,b]\) and let \(x_0, \ldots, x_n\) be a uniform partition of \(I\) with subinterval width \(\Delta x = \frac{b -a}{n}\text{.}\) Consider the subinterval \([x_i, x_{i+1}]\text{.}\) Then

\begin{equation*} \int_{x_i}^{x_{i+1}} f(x) \, dx \approx \frac{f(x_i) + f(x_{i+1})}{2} \Delta x. \end{equation*}

Thus,

\begin{align*} \int_a^b f(x)\, dx \amp = \sum_i \int_{x_i}^{x_{i+1}} f(x) \, dx \\ \amp\approx \sum_i \frac{f(x_i) + f(x_{i+1})}{2} \Delta x \\ \amp = \frac{\Delta x}{2} \left( f(x_0) + 2\sum_{i = 2}^{n-1} f(x_i) + f(x_n)\right)\\ \amp = \frac{b - a}{2n}\left( f(x_0) + 2\sum_{i = 2}^{n-1} f(x_i) + f(x_n)\right) \end{align*}

and we define the \(n\) segement trapezoid approximation to the area under \(f\) by

\begin{equation*} T_n(f, [a,b]) = \frac{b - a}{2n}\left( f(x_0) + 2\sum_{i = 2}^{n-1} f(x_i) + f(x_n)\right). \end{equation*}

Let's use the trapezoid rule to approximate the integral indicated above - \(\int_0^6 x^2 \, dx\text{.}\) We'll use three trapezoids as in the example picture.

We have included the true result, since we can use the fundamental theorem in this case. The result of the approximation with just 3 trapezoids is a relative error of just 5 percent.

Subsection 6.3 A special case of Richardson's extrapolation (optional)

This section will have a different flavor than most of the rest of the notes. Here, we'll see the “analysis” part of numerical analysis - that is, we're going to use theoretical ideas and estimates to improve the approximation given in the trapezoid formula. The idea is that the application of mathematical reasoning can lead to significant improvements in our naive formulations (a theme common in approximation theory).

We'll first recall the triangle inequality, which says that \(\abs{a + b} \leq \abs{a} + \abs{b}\text{.}\) In fact, we can apply this to a sum of any finite length, by induction:

\begin{equation*} \abs{\sum_i a_i} \leq \sum_i \abs{a_i}. \end{equation*}

Now, Riemann sums are finite sums, and so the triangle inequality applies.

\begin{equation*} \abs{\sum_{i = 0}^{n-1} f(x_i^*) \Delta x} \leq \sum_{i = 0}^{n-1} \abs{f(x_i^*)} \Delta x, \end{equation*}

and when the limit of the Riemann sum exists as the number of rectangles tends to infinity (that is, whenever \(f\) is Riemann integrable), we get the integral version of the triangle inequality:

\begin{equation*} \abs{\int_a^b f(x) \, dx} \leq \int_a^b \abs{f(x)} \, dx. \end{equation*}

This is a theorem in real analysis and will be used here without a formal proof beyond the sketch above.

Let \(I = \int_a^b f(x) \, dx\text{.}\) Let \(T_n\) represent the \(n\) segment trapezoid approximation of \(I\text{.}\) Let \(E_T\) be the error in the approximation - that is

\begin{equation*} E_T = I - T_n. \end{equation*}

Our first goal is to measure how large the error \(E_T\) is expected to be in terms of the number of trapezoids \(n\text{.}\)

Let \(\Delta x = x_{i+1} - x_i = \frac{b -a}{n}\text{.}\) We first analyze the error in the trapezoid approximation on a single interval. A \(u\)-substitution gives

\begin{equation*} \int_{x_i}^{x_{i+1}} f(x) \, dx = \int_0^{\Delta x} f(t + x_i) \, dt. \end{equation*}

Using integration by parts twice, we get

\(u\) \(v\)
\(f(t + x_i)\) \(1\)
\(f'(t + x_i)\) \(t + A\)
\(f''(t + x_i)\) \(\frac{(t + A)^2}{2} + B\)
Table 6.3.
where we forgo the usual choice of 0 as the integration constant (hence the \(A, B\) in the table), which gives the formula

\begin{align*} \amp\int_0^{\Delta x} f(t + x_i)\, dx \\ \amp = \left[(t+A)f(t+x_i)\right]_0^{\Delta x} \\ \amp - \left[\left(\frac{(t + A)^2}{2} + B\right)f'(t + x_i)\right]_0^{\Delta x}\\ \amp+ \int_0^{\Delta x} \left(\frac{(t + A)^2}{2} + B\right) f''(t + x_i)\, dx. \end{align*}

From this point, the idea is to choose values of \(A, B\) that force each term to play a certain role, with the goal of concentrating the error of the approximation in the integral term. First, we choose \(A\) so that the first term above is equal to the trapezoid area - that is, we want

\begin{equation*} (\Delta x + A) f(\Delta x + x_i) - Af(x_i) = \frac{f(x_{i+1}) + f(x_i)}{2} \Delta x. \end{equation*}

Algebra shows that \(A = \frac{-\Delta x}{2}\) solves the equation.

Now, we want to choose \(B\) so that the second term is zero. That is, we want

\begin{align*} \amp \left[\left(\frac{(t + A)^2}{2} + B\right)f'(t + x_i)\right]_0^{\Delta x} \\ \amp = \left(\frac{(\Delta x)^2}{8} + B\right)\left[f'(x_{i+1}) - f'(x_i)\right]= 0. \end{align*}

This obviously holds when \(B = \frac{-(\Delta x)^2}{8}\text{.}\)

We conclude that the error on the \(i\)th segment, denoted \(E_T(i)\) is given by

\begin{equation*} E_T(i) = \int_0^{\Delta x} \left(\frac{(t - \frac{\Delta x}{2})^2}{2} - \frac{(\Delta x)^2}{8}\right) f''(t + x_i) \, dt \end{equation*}

Now, we can get the total error in the trapezoid approximation by adding each of the individual errors.

\begin{align*} E_T \amp= \sum E_T(i)\\ \amp = \sum_{i=0}^{n-1} \int_0^{\Delta x} \left(\frac{(t - \frac{\Delta x}{2})^2}{2} - \frac{(\Delta x)^2}{8}\right) f''(t + x_i) \, dt \\ \amp =\int_0^{\Delta x} \left(\frac{(t - \frac{\Delta x}{2})^2}{2} - \frac{(\Delta x)^2}{8}\right) \left(\sum_{i = 0}^{n-1} f''(t + x_i)\right) \, dt \end{align*}

For a well-behaved function \(f\) (the precise assumption is that \(f\) is \(C^2\) on \([a,b]\)),the second derivative is bounded on \([a,b]\) - that is, we assume that there exists a constant \(K\) so that \(\abs{f''(x)} \leq K\) for all \(x \in [a,b]\text{.}\) Then, using the triangle inequality, we derive the approximation

\begin{align*} \abs{E_T} \amp= \abs{\int_0^{\Delta x} \left(\frac{(t - \frac{\Delta x}{2})^2}{2} - \frac{(\Delta x)^2}{8}\right) \left(\sum_{i = 0}^{n-1} f''(t + x_i)\right)\, dt}\\ \amp \leq \int_0^{\Delta x} \abs{\left(\frac{(t - \frac{\Delta x}{2})^2}{2} - \frac{(\Delta x)^2}{8}\right)} \abs{\left(\sum_{i = 0}^{n-1} f''(t + x_i)\right)}\, dt\\ \amp\leq \int_0^{\Delta x} \abs{\left(\frac{(t - \frac{\Delta x}{2})^2}{2} - \frac{(\Delta x)^2}{8}\right)} \abs{\left(\sum_{i = 0}^{n-1} K \right)}\, dt\\ \amp \leq nK \int_0^{\Delta x} \abs{\left(\frac{(t - \frac{\Delta x}{2})^2}{2} - \frac{(\Delta x)^2}{8}\right)} \, dt \end{align*}

The function \(g(t) =\left(\frac{(t - \frac{\Delta x}{2})^2}{2} - \frac{(\Delta x)^2}{8}\right)\) is a parabola that opens upwards with zeros at \(t = 0\) and \(t = h\text{,}\) and so

\begin{align*} \amp \int_0^{\Delta x} \abs{\left(\frac{(t - \frac{\Delta x}{2})^2}{2} - \frac{(\Delta x)^2}{8}\right)} \, dt\\ \amp \int_0^{\Delta x} \frac{(\Delta x)^2}{8} - \left(\frac{(t - \frac{\Delta x}{2})^2}{2} \right) \, dt\\ \amp = \left[\frac{(\Delta x)^2}{8}t - \frac{(t - \frac{\Delta x}{2})^3}{6}\right]_0^{\Delta x} \\ \amp= \frac{(\Delta x)^3}{12} \end{align*}

Putting this together with the previous computation, we get

\begin{equation*} \abs{E_T} \leq nK\frac{(\Delta x)^3}{12} = nK \frac{(b-a)^3}{12 n^3} = \frac{K(b-a)^3}{12 n^2}, \end{equation*}

where \(K\) was an absolute bound for \(f''\) on \([a,b]\text{,}\) and in the worst case, we get \(\abs{E_T} \leq \frac{K(b-a)^3}{12} \frac{1}{n^2} = \frac{C}{n^2}\) - that is, the error is proportional to \(\frac{1}{n^2}\text{,}\) which establishes the claim.

Under the assumption of worst case error and a reasonable function \(f\text{,}\) we conclude that the total trapezoidal error \(E_T\) is proportional to \(\frac{1}{n^2}\text{,}\) or in other words that

\begin{equation*} E_T = \frac{C}{n^2}. \end{equation*}

So how can we use this to build a better process? Note that for \(n\) segments, we can write

\begin{equation*} I = T_n + \frac{C}{n^2} \end{equation*}

and likewise for \(2n\) segments, we have

\begin{equation*} I = T_{2n} + \frac{C}{(2n)^2}, \end{equation*}

which is a system of simultaneous equations. We'll prepare to eliminate \(C\text{.}\)

\begin{align*} I \amp = T_{n} + \frac{C}{n^2}\\ I \amp = T_{2n} + \frac{C}{4n^2}\\ \amp\\ n^2 I \amp = n^2 T_n + C \\ 4n^2 I \amp= 4n^2 T_{2n} + C \\ \amp \\ 3n^2 I \amp = 4n^2 T_{2n} - n^2 T_n\\ \amp \\ I \amp = T_{2n} + \frac{T_{2n} - T_n}{3} \end{align*}

Thus, we have what is known as a first order Richardson's extrapolation -

\begin{equation*} \int_a^b f(x) \, dx \approx T_{2n} + \frac{T_{2n} - T_n}{3}. \end{equation*}

Let's see how it performs with our existing example.

Subsection 6.4 Simpson's 1/3 rule

An alternative to using trapezoids is to use polynomials to interpolate sample points. It turns out that using quadratic polynomials on equally spaced interpolation points gives a very nice formula. We'll begin with a single segment and approximate \(\int_a^b f(x)\, dx\text{.}\)

Recall that there is a unique parabola through any three points - we'll use the points \((a, f(a)), ((a+b)/2, f((a+b)/2), (b, f(b))\text{.}\) We have several techniques available for finding such an interpolation - we'll derive ours using Newton polynomials. The Newton polynomial through our points is

\begin{equation*} g(x) = b_0 + b_1 (x - a) + b_2 (x -a) (x - (\frac{a+b}{2})), \end{equation*}

where

\begin{equation*} b_0 = a, b_1 = \frac{f((a+b)/2) - f(a)}{(a + b)/2 - a}, b_2 = \frac{\frac{f(b) - f((a+b)/2)}{b - (a+b)/2} - \frac{f((a+b)/2) - f(a)}{(a+b)/2 - a}}{b - a}. \end{equation*}

So in effect, we're saying that

\begin{align*} \int_a^b f(x)\, dx \amp \approx \int_a^b b_0 + b_1 (x - a) + b_2 (x -a) (x - (\frac{a+b}{2}))\\ \amp = \ldots\\ \amp= \frac{\Delta x}{3}\left( f(a) + 4f(\frac{a+b}{2}) + f(b)\right) \end{align*}

where we have shoved integration and substitution under the rug (or into the ellipsis, as it were).

Now suppose that \([a,b]\) is partitioned into \(n\) segments of equal length \(\Delta x = \frac{b - a}{n}\text{,}\) where \(n\) is even (which will allow overlapping sequences of three points). Then for three sequential points \(x_i, x_{i+1}, x_{i+2}\text{,}\) the previous computation gives

\begin{equation*} \int_{x_i}^{x_{i+2}} f(x) \, dx \approx \frac{\Delta x}{3}\left(f(x_i) + 4f(x_{i+1}) + f(x_{i+2})\right). \end{equation*}

Since we chose \(n\) to be even, we get

\begin{align*} \int_a^b f(x)\, dx \amp\approx \frac{\Delta x}{3}\left(f(x_0) + 4f(x_{1}) + f(x_{2})\right) + \ldots + \frac{\Delta x}{3}\left(f(x_{n-2}) + 4f(x_{n-1}) + f(x_{n})\right)\\ \amp\approx \frac{\Delta x}{3} \left(f(x_0) + 4 \sum_{i \text{ odd}} f(x_i)) + 2 \sum_{i \text{ even}} f(x_i) + f(x_n)\right) \end{align*}
login