Numerical Integration

Earlier in this text we defined the definite integral of a function over an interval as the limit of Riemann sums. In general, any Riemann sum of a function \(f(x)\) over an interval \([a,b]\) may be viewed as an estimate of \({\int }_{a}^{b}f(x)dx.\) Recall that a Riemann sum of a function \(f(x)\) over an interval \([a,b]\) is obtained by selecting a partition

and a set

The Riemann sum corresponding to the partition \(P\) and the set \(S\) is given by \(\sum _{i=1}^{n}f({x}_{i}^{*})\text{\Delta }{x}_{i},\) where \(\text{\Delta }{x}_{i}={x}_{i}-{x}_{i-1},\) the length of the ith subinterval.

The midpoint rule for estimating a definite integral uses a Riemann sum with subintervals of equal width and the midpoints, \({m}_{i},\) of each subinterval in place of \({x}_{i}^{*}.\) Formally, we state a theorem regarding the convergence of the midpoint rule as follows.

As we can see in , if \(f(x)\ge 0\) over \([a,b],\) then \(\sum _{i=1}^{n}f({m}_{i})\text{\Delta }x\) corresponds to the sum of the areas of rectangles approximating the area between the graph of \(f(x)\) and the x-axis over \([a,b].\) The graph shows the rectangles corresponding to \({M}_{4}\) for a nonnegative function over a closed interval \([a,b].\)

Condensed — the full section is in OpenStax Calculus Volume 2.

We can also approximate the value of a definite integral by using trapezoids rather than rectangles. In , the area beneath the curve is approximated by trapezoids rather than by rectangles.

The trapezoidal rule for estimating definite integrals uses trapezoids rather than rectangles to approximate the area under a curve. To gain insight into the final form of the rule, consider the trapezoids shown in . We assume that the length of each subinterval is given by \(\text{\Delta }x.\) First, recall that the area of a trapezoid with a height of h and bases of length \({b}_{1}\) and \({b}_{2}\) is given by \(\text{Area}=\frac{1}{2}h({b}_{1}+{b}_{2}).\) We see that the first trapezoid has a height \(\text{\Delta }x\) and parallel bases of length \(f({x}_{0})\) and \(f({x}_{1}).\) Thus, the area of the first trapezoid in is

The areas of the remaining three trapezoids are

Consequently,

After taking out a common factor of \(\frac{1}{2}\text{\Delta }x\) and combining like terms, we have

Generalizing, we formally state the following rule.

Before continuing, let’s make a few observations about the trapezoidal rule. First of all, it is useful to note that

Condensed — the full section is in OpenStax Calculus Volume 2.

An important aspect of using these numerical approximation rules consists of calculating the error in using them for estimating the value of a definite integral. We first need to define absolute error and relative error.

In the two previous examples, we were able to compare our estimate of an integral with the actual value of the integral; however, we do not typically have this luxury. In general, if we are approximating an integral, we are doing so because we cannot compute the exact value of the integral itself easily. Therefore, it is often helpful to be able to determine an upper bound for the error in an approximation of an integral. The following theorem provides error bounds for the midpoint and trapezoidal rules. The theorem is stated without proof.

We can use these bounds to determine the value of \(n\) necessary to guarantee that the error in an estimate is less than a specified value.

Condensed — the full section is in OpenStax Calculus Volume 2.

With the midpoint rule, we estimated areas of regions under curves by using rectangles. In a sense, we approximated the curve with piecewise constant functions. With the trapezoidal rule, we approximated the curve by using piecewise linear functions. What if we were, instead, to approximate a curve using piecewise quadratic functions? With Simpson’s rule, we do just this. We partition the interval into an even number of subintervals, each of equal width. Over the first pair of subintervals we approximate \({\int }_{{x}_{0}}^{{x}_{2}}f(x)dx\) with \({\int }_{{x}_{0}}^{{x}_{2}}p(x)dx,\) where \(p(x)=A{x}^{2}+Bx+C\) is the quadratic function passing through \(({x}_{0},f({x}_{0})),\) \(({x}_{1},f({x}_{1})),\) and \(({x}_{2},f({x}_{2}))\) (). Over the next pair of subintervals we approximate \({\int }_{{x}_{2}}^{{x}_{4}}f(x)dx\) with the integral of another quadratic function passing through \(({x}_{2},f({x}_{2})),\) \(({x}_{3},f({x}_{3})),\) and \(({x}_{4},f({x}_{4})).\) This process is continued with each successive pair of subintervals.

To understand the formula that we obtain for Simpson’s rule, we begin by deriving a formula for this approximation over the first two subintervals. As we go through the derivation, we need to keep in mind the following relationships:

\({x}_{2}-{x}_{0}=2\text{\Delta }x,\) where \(\text{\Delta }x\) is the length of a subinterval.

Thus,

If we approximate \({\int }_{{x}_{2}}^{{x}_{4}}f(x)dx\) using the same method, we see that we have

Combining these two approximations, we get

The pattern continues as we add pairs of subintervals to our approximation. The general rule may be stated as follows.

Condensed — the full section is in OpenStax Calculus Volume 2.

• We can use numerical integration to estimate the values of definite integrals when a closed form of the integral is difficult to find or when an approximate value only of the definite integral is needed.
• The most commonly used techniques for numerical integration are the midpoint rule, trapezoidal rule, and Simpson’s rule.
• The midpoint rule approximates the definite integral using rectangular regions whereas the trapezoidal rule approximates the definite integral using trapezoidal approximations.
• Simpson’s rule approximates the definite integral by first approximating the original function using piecewise quadratic functions.

Approximate the following integrals using either the midpoint rule, trapezoidal rule, or Simpson’s rule as indicated. (Round answers to three decimal places.)

Approximate the integral to three decimal places using the indicated rule.

When we first explored finding the net signed area bounded by a curve, we developed the concept of a Riemann sum as a helpful estimation tool and a key step in the definition of the definite integral. Recall that the left, right, and middle Riemann sums of a function \(f\) on an interval \([a,b]\) are given by \[\begin{aligned}L_n = f(x_0) \Delta x + f(x_1) \Delta x + \cdots + f(x_{n-1}) \Delta x \amp= \sum_{i = 0}^{n-1} f(x_i) \Delta x, \\ R_n = f(x_1) \Delta x + f(x_2) \Delta x + \cdots + f(x_{n}) \Delta x \amp= \sum_{i = 1}^{n} f(x_i) \Delta x, \\ M_n = f(\overline{x}_1) \Delta x + f(\overline{x}_2) \Delta x + \cdots + f(\overline{x}_{n}) \Delta x \amp= \sum_{i = 1}^{n} f(\overline{x}_i) \Delta x\end{aligned}\], where \(x_0 = a\), \(x_i = a + i\Delta x\), \(x_n = b\), and \(\Delta x = \frac{b-a}{n}\). For the middle sum, we defined \(\overline{x}_{i} = (x_{i-1} + x_i)/2\).

A Riemann sum is a sum of (possibly signed) areas of rectangles. The value of \(n\) determines the number of rectangles, and our choice of left endpoints, right endpoints, or midpoints determines the heights of the rectangles. We can see the similarities and differences among these three options in Figure, where we consider the function \(f(x) = \frac{1}{20}(x-4)^3 + 7\) on the interval \([1,8]\), and use 5 rectangles for each of the Riemann sums.

While it is a good exercise to compute a few Riemann sums by hand, just to ensure that we understand how they work and how varying the function, the number of subintervals, and the choice of endpoints or midpoints affects the result, using computing technology is the best way to determine \(L_n\), \(R_n\), and \(M_n\). Any computer algebra system will offer this capability; as we saw in Preview Activity, a straightforward option that is freely available online is this interactive. Note that we can adjust the formula for \(f(x)\), the window of \(x\)- and \(y\)-values of interest, the number of subintervals, and the type of Riemann sum. (See Preview Activity for any needed reminders on how the interactive works.)

In this section we explore several different alternatives for estimating definite integrals. Our main goal is to develop formulas to estimate definite integrals accurately without using a large numbers of rectangles.

So far, we have used the simplest possible quadrilaterals (that is, rectangles) to estimate areas. It is natural, however, to wonder if other familiar shapes might serve us even better.

An alternative to \(L_n\), \(R_n\), and \(M_n\) is called the Trapezoid Rule. Rather than using a rectangle to estimate the (signed) area bounded by \(y = f(x)\) on a small interval, we use a trapezoid. For example, in Figure, we estimate the area under the curve using three subintervals and the trapezoids that result from connecting the corresponding points on the curve with straight lines.

The biggest difference between the Trapezoid Rule and a Riemann sum is that on each subinterval, the Trapezoid Rule uses two function values, rather than one, to estimate the (signed) area bounded by the curve. For instance, to compute \(D_1\), the area of the trapezoid on \([x_0, x_1]\), we observe that the left base has length \(f(x_0)\), while the right base has length \(f(x_1)\). The height of the trapezoid is \(x_1 - x_0 = \Delta x = \frac{b-a}{3}\). The area of a trapezoid is the average of the bases times the height, so we have \[\begin{aligned}\end{aligned}\].

Using similar computations for \(D_2\) and \(D_3\), we find that \(T_3\), the trapezoidal approximation to \(\int_a^b f(x) \, dx\) is given by \[\begin{aligned}T_3 \amp= D_1 + D_2 + D_3 \\ \amp= \frac{1}{2}(f(x_0) + f(x_1)) \cdot \Delta x + \frac{1}{2}(f(x_1) + f(x_2)) \cdot \Delta x + \frac{1}{2}(f(x_2) + f(x_3)) \cdot \Delta x\end{aligned}\].

Because both left and right endpoints are being used, we recognize within the trapezoidal approximation the use of both left and right Riemann sums. Rearranging the expression for \(T_3\) by removing factors of \(\frac{1}{2}\) and \(\Delta x\), grouping the left endpoint and right endpoint evaluations of \(f\), we see that \[\begin{aligned}\end{aligned}\].

We now observe that two familiar sums have arisen. The left Riemann sum \(L_3\) is \(L_3 = f(x_0) \Delta x + f(x_1) \Delta x + f(x_2) \Delta x\), and the right Riemann sum is \(R_3 = f(x_1) \Delta x + f(x_2) \Delta x + f(x_3) \Delta x\). Substituting \(L_3\) and \(R_3\) for the corresponding expressions in Equation, it follows that \(T_3 = \frac{1}{2} \left[ L_3 + R_3 \right]\). We have thus seen a very important result: using trapezoids to estimate the (signed) area bounded by a curve is the same as averaging the estimates generated by using left and right endpoints.

Condensed — the full section is in Boelkins, Active Calculus.

We know from the definition of the definite integral that if we let \(n\) be large enough, we can make any of the approximations \(L_n\), \(R_n\), and \(M_n\) as close as we'd like (in theory) to the exact value of \(\int_a^b f(x) \, dx\). Thus, it may be natural to wonder why we ever use any rule other than \(L_n\) or \(R_n\) (with a sufficiently large \(n\) value) to estimate a definite integral. One of the primary reasons is that as \(n \to \infty\), \(\Delta x = \frac{b-a}{n} \to 0\), and thus in a Riemann sum calculation with a large \(n\) value, we end up multiplying by a number that is very close to zero. Doing so often generates roundoff error, because representing numbers close to zero accurately is a persistent challenge for computers.

Hence, we explore ways to estimate definite integrals to high levels of precision, but without using extremely large values of \(n\). Paying close attention to patterns in errors, such as those observed in Activity, is one way to begin to see some alternate approaches.

To begin, we compare the errors in the Midpoint and Trapezoid rules. First, consider a function that is concave up on a given interval, and picture approximating the area bounded on that interval by both the Midpoint and Trapezoid rules using a single subinterval, as shown in Figure.

As seen at left in Figure, it is evident that whenever the function is concave up on an interval, the Trapezoid Rule with one subinterval, \(T_1\), will overestimate the exact value of the definite integral on that interval. From a careful analysis of the line that bounds the top of the rectangle for the Midpoint Rule (shown in magenta), we see that if we rotate this line segment until it is tangent to the curve at the midpoint of the interval (as shown at right in Figure), the resulting trapezoid has the same area as \(M_1\), and this value is less than the exact value of the definite integral. Thus, when the function is concave up on the interval, \(M_1\) underestimates the integral's true value.

These observations extend easily to the situation where the function's concavity remains consistent but we use larger values of \(n\) in the Midpoint and Trapezoid Rules. Hence, whenever \(f\) is concave up on \([a,b]\), \(T_n\) will overestimate the value of \(\int_a^b f(x) \, dx\), while \(M_n\) will underestimate \(\int_a^b f(x) \, dx\). The reverse observations are true in the situation where \(f\) is concave down.

Condensed — the full section is in Boelkins, Active Calculus.

When we first developed the Trapezoid Rule, we observed that it is an average of the Left and Right Riemann sums: \[\begin{aligned}\end{aligned}\].

If a function is always increasing or always decreasing on the interval \([a,b]\), one of \(L_n\) and \(R_n\) will over-estimate the true value of \(\int_a^b f(x) \, dx\), while the other will under-estimate the integral. Thus, the errors that result from \(L_n\) and \(R_n\) will have opposite signs; so averaging \(L_n\) and \(R_n\) eliminates a considerable amount of the error present in the respective approximations. In a similar way, it makes sense to think about averaging \(M_n\) and \(T_n\) in order to generate a still more accurate approximation.

We've just observed that \(M_n\) is typically about twice as accurate as \(T_n\). So we use the weighted average \[\begin{aligned}\end{aligned}\].

The rule for \(S_{2n}\) giving by Equation is usually known as Simpson's Rule. Thomas Simpson was an 18th century mathematician; his idea was to extend the Trapezoid rule, but rather than using straight lines to build trapezoids, to use quadratic functions to build regions whose area was bounded by parabolas (whose areas he could find exactly). Simpson's Rule is often developed from the more sophisticated perspective of using interpolation by quadratic functions. Note that we use \(S_{2n}\) rather that \(S_n\) since the \(n\) points the Midpoint Rule uses are different from the \(n\) points the Trapezoid Rule uses, and thus Simpson's Rule is using \(2n\) points at which to evaluate the function. We build upon the results in Table to see the approximations generated by Simpson's Rule. In particular, in Table, we include all of the results in Table, but include additional results for \(S_8 = \frac{2M_4 + T_4}{3}\) and \(S_{16} = \frac{2M_8 + T_8}{3}\).

Condensed — the full section is in Boelkins, Active Calculus.

As we conclude our discussion of numerical approximation of definite integrals, it is important to summarize general trends in how the various rules over- or under-estimate the true value of a definite integral, and by how much. To revisit some past observations and see some new ones, we consider the following activity.

The results seen in Activity generalize nicely. For instance, if \(f\) is decreasing on \([a,b]\), \(L_n\) will over-estimate the exact value of \(\int_a^b f(x) \,dx\), and if \(f\) is concave down on \([a,b]\), \(M_n\) will over-estimate the exact value of the integral. An excellent exercise is to write a collection of scenarios of possible function behavior, and then categorize whether each of \(L_n\), \(R_n\), \(T_n\), and \(M_n\) is an over- or under-estimate.

Finally, we make two important notes about Simpson's Rule. When T.Simpson first developed this rule, his idea was to replace the function \(f\) on a given interval with a quadratic function that shared three values with the function \(f\). In so doing, he guaranteed that this new approximation rule would be exact for the definite integral of any quadratic polynomial. In one of the pleasant surprises of numerical analysis, it turns out that even though it was designed to be exact for quadratic polynomials, Simpson's Rule is exact for any cubic polynomial: that is, if we are interested in an integral such as \(\int_2^5 (5x^3 - 2x^2 + 7x - 4)\, dx\), \(S_{2n}\) will always be exact, regardless of the value of \(n\). This is just one more piece of evidence that shows how effective Simpson's Rule is as an approximation tool for estimating definite integrals. One reason that Simpson's Rule is so effective is that \(S_{2n}\) benefits from using \(2n+1\) points of data. Because it combines \(M_n\), which uses \(n\) midpoints, and \(T_n\), which uses the \(n+1\) endpoints of the chosen subintervals, \(S_{2n}\) takes advantage of the maximum amount of information we have when we know function values at the endpoints and midpoints of \(n\) subintervals.

• For a definite integral such as \(\int_0^1 e^{-x^2} \, dx\) when we cannot use the First Fundamental Theorem of Calculus because the integrand lacks an elementary algebraic antiderivative, we can estimate the integral's value by using a sequence of Riemann sum approximations. Typically, we start by computing \(L_n\), \(R_n\), and \(M_n\) for one or more chosen values of \(n\).

• The Trapezoid Rule, which estimates \(\int_a^b f(x) \, dx\) by using trapezoids, rather than rectangles, can also be viewed as the average of Left and Right Riemann sums. That is, \(T_n = \frac{1}{2}(L_n + R_n)\).

• The Midpoint Rule is typically twice as accurate as the Trapezoid Rule, and the signs of the respective errors of these rules are opposites. Hence, by taking the weighted average \(S_{2n} = \frac{2M_n + T_n}{3}\), we can build a much more accurate approximation to \(\int_a^b f(x) \, dx\) by using approximations we have already computed. The rule for \(S_{2n}\) is known as Simpson's Rule, which can also be developed by approximating a given continuous function with pieces of quadratic polynomials.