Numerical Integration

We have seen that definite integrals arise in many different areas and that the Fundamental Theorem of Calculus is a powerful tool for evaluating definite integrals. However, it cannot always be applied: there are some functions which do not have an antiderivative which can be expressed in terms of familiar functions such as polynomials, exponentials and trigonometric functions. One such example is $ e^{-x^2} $ ; of course, this is an important function since it is the probability density function for the normal distribution.

In this section, we will demonstrate two tools for approximating a definite integral to any degree of accuracy. These tools can be used to evaluate definite integrals when the integrand has no simple antiderivative. Moreover, we sometimes only have information about a function by making observations at a certain number of points. In that case, we do not have a nice formula for the function we are integrating, but only some data points. Our methods can help us to evaluate a definite integral in this case as well.


The Trapezoidal Rule

In our first method, we will approximate the definite integral $ \int_a^b f(x) ~dx $ by approximating the graph of the function by a straight line. That means that we approximate the area under the graph $ y = f(x) $ by the trapezoid formed below.

The area of the trapezoid is just the area of the rectangle plus the area of the triangle. That means our approximation is

\begin{eqnarray*} \int_a^b f(x)~dx & = & (b-a)f(a) + \frac 12 (b-a)\big[f(b) - f(a)\big] \\ & = & (b-a)\big[ f(a) + \frac 12 f(b) - \frac 12 f(a)\big] \\ & = & \frac 12(b-a) (f(a) + f(a)) \end{eqnarray*}

Of course, we cannot expect this to be a good approximation. However, we can break the region $ [a,b] $ into many smaller pieces and apply the approximation on each piece. On the smaller pieces, the graph looks more and more like a straight line so the approximation should improve.

Let's choose some positive integer n and break the interval $ [a,b] $ into n equal pieces. The width of each piece is $ h = \frac{b-a}{n} $ .

We will label the points defined by the sub-intervals by $ x_i $ and call $ y_i = f(x_i) $ . If we the approximate the area under the graph by the area of the trapezoids, we have

\begin{eqnarray*} \int_a^b f(x)~dx \approx T_n & = & \frac h2(y_0 + y_1) + \frac h2(y_1 + y_2) + \ldots | \frac h2(y_{n-1} + y_n) \\ & = & \frac h2\big[y_0 + y_1 + y_1 + y_2 + y_2 + y_3 + \ldots + y_n \big] \\ & = & \frac h2 \big[ y_0 + 2y_1 + 2y_2 + \ldots + 2y_{n-1} + y_n\big] \end{eqnarray*}

Example: We will consider the integral $ \int_0^1 (1-x^2) ~dx $ . Of course, we can evaluate this integral directly with the Fundamental Theorem of Calculus: $ \int_0^1(1-x^2)~dx = (x-\frac{x^3}{3})|_0^1 = \frac 23 $ . We now will build the trapezoidal approximation so that we might see explicitly how much of an error we are making.

First, we will consider $ n = 1 $ . In this case, $ h = \frac {1-0}{1} = 1 $ . There are only two points $ x_0 = 0 $ and $ x_1 = 1 $ . Then we have $ y_0 = 1-x_0^2 = 1 $ and $ y_1 = 1-x_1^2 = 0 $ . Then the trapezoidal rule produces

\[ \int_0^1 (1-x^2)~dx \approx T_1 = \frac h2 (y_0 + y_1) = \frac 12(1+0) =\frac 12 \]

This means that the Trapezoidal Rule with $ n = 1 $ produces an approximation of $ \frac 12 $ to the integral which we know is $ \frac 23 $ .

Now let's see what happens when $ n = 2 $ . In this case, we have $ h = \frac{1-0}{2} = \frac 12 $ . The points for us to consider are $ x_0 = 0, x_1 = \frac 12 $ and $ x_2 =1 $ . This produces $ y_0 = 1, y_1 = \frac 34 $ and $ y_2 = 0 $ . Then we have

\[ \int_0^1(1-x^2)~dx \approx T_2 = \frac h2(y_0 + 2y_1 + y_2) = \frac 14(1+2\frac 34 + 0) = \frac 58 = 0.625. \]

Clearly, this is a better approximation. The following demonstration will show what happens when we increase n further. Notice that the approximation becomes quite good. The reason for this is clear from the picture: as the size of the intervals becomes very small, the graph is better approximated by a straight line on each interval.


Simpson's Rule

In the example above, you can see that the Trapezoidal Rule provides a reasonable approximation to a definite integral if we take a large number of steps. Notice that the error in the approximation originates in the fact that general graphs are curved and we are approximating them by straight lines. We will now form an approximation which takes into account the curvature of the graph: the result is a more efficient approximation called Simpson's Rule.

Simpson's Rule is formed by approximating a general curve by a parabola.

In this picture, the red graph is a parabola which approximates the blue graph. Remember that a parabola is the graph of a quadratic function $ y = ax^2 + bx + c $ and so three pieces of information are required to determine the coefficients $ a,b $ , and $ c $ . This means that we must use three data points $ (x_0,y_0), (x_1,y_1) $ and $ (x_2,y_2) $ to fix the parabola.

We won't show you how to determine the coefficients for the parabola, but it is fairly straightforward. As before, the width of each of the two intervals is $ h = \frac{b-a}{2} $ . With a little bit of work, you would find the approximation

\[ \int_a^b f(x)~dx \approx S_1 = \frac h3\big[y_0 + 4y_1 + y_2\big] \]

This is Simpson's Rule with one step. More generally, we can break the interval into several pieces and apply Simpson's Rule on each interval. For instance, to use n steps, break the interval $ [a,b] $ into $ 2n $ pieces, each of width $ h = \frac{b-a}{2n} $ . Call the x coordinates $ x_0,x_1,\ldots,x_{2n-1},x_{2n} $ and let $ y_i=f(x_i) $ . Then we have

\[ \int_a^bf(x)~dx \approx S_n = \frac h3\big[y_0 + 4y_1 + 2y_2 + 4y_3 + 2y_4 +\ldots + 4y_{2n-1} + y_{2n}\big] \]

This is the same idea as the Trapezoidal Rule but you see that the algorithm is slightly different. Inside the sum, the endpoints are weighted once, while the odd values of $ y $ are weigted four times and the even values of $ y $ in the middle are weighted twice.

Examples:

  1. First, we'll consider the same example as we studied with the Trapezoidal Rule: $ \int_0^1 (1-x^2)~dx = \frac 23 $ .

    Let's apply Simpson's Rule with one term: $ n = 1 $ so that $ h = \frac 12 $ . Then $ x_0 = 1, x_1 = \frac 12 $ and $ x_2 = 1 $ . This means that $ y_0 = 1, y_1 = \frac 34 $ and $ y_2 = 0 $ .

    Now applying Simpson's Rule gives

    \[ \int_0^1 (1-x^2)~dx \approx S_1 = \frac 16\big[1 + 4\frac 34 + 0\big] = \frac 23 \]

    This means that, with one step, Simpson's Rule gives the correct answer. This shouldn't be too surprising since the Rule uses a parabola to approximate the graph. In this case, the graph we are interested in is already a parabola and so Simpson's Rule will produce the correct answer. In comparing to the Trapezoidal Rule, you can see the savings in effort: there we needed many steps to get a good approximation to the integral.

  2. Consider the integral $ \int_0^1 \frac{1}{1+x^2}~dx = \tan^{-1} x|_0^1 = \frac \pi 4 $ .

    You have probably been told for most of your mathematical career that $ \pi \approx 3.14159... $ but you have probably never seen where that number comes from. We can use Simpson's Rule to approximate the integral and then use the fact that

    \[ \pi = 4\int_0^1\frac{1}{1+x^2}~dx \]

    to obtain an approximation for $ \pi $ .

    With $ n = 1 $ , we have $ h = \frac 12 $ and $ x_0 = 0,x_1 = \frac 12 $ and $ x_2 = 1 $ . This means that $ y_0 = 1,y_1 = \frac 45 $ and $ y_2 = \frac 12 $ .

    Our approximation is then

    \[ \int_0^1 \frac{1}{1+x^2}~dx \approx S_1 = \frac 16(1 + 4\frac 45 + \frac 12) = \frac {47}{60} \]

    This gives an approximation for $ \pi \approx 4\frac {47}{60} = \frac {47}{15} = 3.1333... $ . So with just one step, we have a pretty good approximation for $ \pi $ .

    With $ n = 2 $ , we have

    \[ \int_0^1 \frac{1}{1+x^2}~dx \approx S_2 = \frac{1}{12}\big[1 + 4\frac{1}{1+(\frac 14)^2} + 2 \frac{1}{1+(\frac12)^2} + 4\frac{1}{1+(\frac 34)^2} + \frac{1}{1+1^2}\big] = 0.7853 \]

    This produces an approximation $ \pi \approx 4\int_0^1\frac{1}{1+x^2}~dx =3.141568. $ With just two steps of Simpson's Rule, we have computed $ \pi $ to four decimal places!

    The demonstration below will provide approximations for this integral using both the Trapezoidal Rule and Simpson's Rule. You can see how much better Simpson's Rule really is.

  3. Of course, the real value of Simpson's Rule is in computing integrals of functions which cannot be antidifferentiated using familiar functions. One such example is

    \[ \int_{-1}^1 e^{-x^2}~dx \approx 1.4936 \]

    where the approximation is obtained using Simpson's Rule. This integral is needed to determine the probability that a typical sample lies within one standard deviation of the mean in the normal distribution.