Thoughts on Numerical Integration (Part 16): Midpoint rule and local rate of convergence

Numerical integration is a standard topic in first-semester calculus. From time to time, I have received questions from students on various aspects of this topic, including:

Why is numerical integration necessary in the first place?
Where do these formulas come from (especially Simpson’s Rule)?
How can I do all of these formulas quickly?
Is there a reason why the Midpoint Rule is better than the Trapezoid Rule?
Is there a reason why both the Midpoint Rule and the Trapezoid Rule converge quadratically?
Is there a reason why Simpson’s Rule converges like the fourth power of the number of subintervals?

In this series, I hope to answer these questions. While these are standard questions in a introductory college course in numerical analysis, and full and rigorous proofs can be found on Wikipedia and Mathworld, I will approach these questions from the point of view of a bright student who is currently enrolled in calculus and hasn’t yet taken real analysis or numerical analysis.

In this post, we will perform an error analysis for the Midpoint Rule

$\int_a^b f(x) \, dx \approx h \left[f(c_1) + f(c_2) + \dots + f(c_n) \right] \equiv M_n$

where $n$ is the number of subintervals and $h = (b-a)/n$ is the width of each subinterval, so that $x_k = x_0 + kh$ . Also, $c_i = (x_{i-1} + x_i)/2$ is the midpoint of the $i$ th subinterval.

As noted above, a true exploration of error analysis requires the generalized mean-value theorem, which perhaps a bit much for a talented high school student learning about this technique for the first time. That said, the ideas behind the proof are accessible to high school students, using only ideas from the secondary curriculum (especially the Binomial Theorem), if we restrict our attention to the special case $f(x) = x^k$ , where $k \ge 5$ is a positive integer.

For this special case, the true area under the curve $f(x) = x^k$ on the subinterval $[x_i, x_i +h]$ will be

$\displaystyle \int_{x_i}^{x_i+h} x^k \, dx = \frac{1}{k+1} \left[ (x_i+h)^{k+1} - x_i^{k+1} \right]$

$= \displaystyle \frac{1}{k+1} \left[x_i^{k+1} + {k+1 \choose 1} x_i^k h + {k+1 \choose 2} x_i^{k-1} h^2 + {k+1 \choose 3} x_i^{k-2} h^3 + {k+1 \choose 4} x_i^{k-3} h^4+ {k+1 \choose 5} x_i^{k-4} h^5+ O(h^6) - x_i^{k+1} \right]$

$= \displaystyle \frac{1}{k+1} \bigg[ (k+1) x_i^k h + \frac{(k+1)k}{2} x_i^{k-1} h^2 + \frac{(k+1)k(k-1)}{6} x_i^{k-2} h^3+ \frac{(k+1)k(k-1)(k-2)}{24} x_i^{k-3} h^4$

$+ \displaystyle \frac{(k+1)k(k-1)(k-2)(k-3)}{120} x_i^{k-4} h^5 \bigg] + O(h^6)$

$= x_i^k h + \displaystyle \frac{k}{2} x_i^{k-1} h^2 + \frac{k(k-1)}{6} x_i^{k-2} h^3 + \frac{k(k-1)(k-2)}{24} x_i^{k-3} h^4 + \frac{k(k-1)(k-2)(k-3)}{120} x_i^{k-4} h^5 + O(h^6)$

In the above, the shorthand $O(h^6)$ can be formally defined, but here we’ll just take it to mean “terms that have a factor of $h^6$ or higher that we’re too lazy to write out.” Since $h$ is supposed to be a small number, these terms will small in magnitude and thus can be safely ignored. I wrote the above formula to include terms up to and including $h^5$ because I’ll need this later in this series of posts. For now, looking only at the Midpoint Rule, it will suffice to write this integral as

$\displaystyle \int_{x_i}^{x_i+h} x^k \, dx =x_i^k h + \displaystyle \frac{k}{2} x_i^{k-1} h^2 + \frac{k(k-1)}{6} x_i^{k-2} h^3 + O(h^4)$ .

Using the midpoint of the subinterval, the left-endpoint approximation of $\displaystyle \int_{x_i}^{x_i+h} x^k \, dx$ is $\displaystyle \left(x_i+ \frac{h}{2} \right)^k h$ . Using the Binomial Theorem, this expands as

$x_i^k h + \displaystyle {k \choose 1} x_i^{k-1} \frac{h^2}{2} + {k \choose 2} x_i^{k-2} \frac{h^3}{4} + {k \choose 3} x_i^{k-3} \frac{h^4}{8} + {k \choose 4} x_i^{k-4} \frac{h^5}{16} + O(h^6)$

$= x_i^k h + \displaystyle \frac{k}{2} x_i^{k-1} h^2 + \frac{k(k-1)}{8} x_i^{k-2} h^3 + \frac{k(k-1)(k-2)}{48} x_i^{k-3} h^4$

$\displaystyle + \frac{k(k-1)(k-2)(k-3)}{384} x_i^{k-4} h^5 + O(h^6)$

Once again, this is a little bit overkill for the present purposes, but we’ll need this formula later in this series of posts. Truncating somewhat earlier, we find that the Midpoint Rule for this subinterval gives

$x_i^k h + \displaystyle \frac{k}{2} x_i^{k-1} h^2 + \displaystyle \frac{k(k-1)}{8} x_i^{k-2} h^3 + O(h^4)$

Subtracting from the actual integral, the error in this approximation will be equal to

$\displaystyle x_i^k h + \frac{k}{2} x_i^{k-1} h^2 + \frac{k(k-1)}{6} x_i^{k-2} h^3 - x_i^k h - \frac{k}{2} x_i^{k-1} h^2 - \frac{k(k-1)}{8} x_i^{k-2} h^3 + O(h^4)$

$= \displaystyle \frac{k(k-1)}{24} x_i^{k-2} h^3 + O(h^4)$

In other words, unlike the left-endpoint and right-endpoint approximations, both of the first two terms $x_i^k h$ and $\displaystyle \frac{k}{2} x_i^{k-1} h^2$ cancel perfectly, leaving us with a local error on the order of $h^3$ .

The logic for determining the global error is much the same as what we used earlier for the left-endpoint rule. The total error when approximating $\displaystyle \int_a^b x^k \, dx = \int_{x_0}^{x_n} x^k \, dx$ will be the sum of the errors for the integrals over $[x_0,x_1]$ , $[x_1,x_2]$ , through $[x_{n-1},x_n]$ . Therefore, the total error will be

$E \approx \displaystyle \frac{k(k-1)}{24} \left(x_0^{k-2} + x_1^{k-2} + \dots + x_{n-1}^{k-2} \right) h^3$ .

So that this formula doesn’t appear completely mystical, this actually matches the numerical observations that we made earlier. The figure below shows the left-endpoint approximations to $\displaystyle \int_1^2 x^9 \, dx$ for different numbers of subintervals. If we take $n = 100$ and $h = 0.01$ , then the error should be approximately equal to

$\displaystyle \frac{9 \times 8}{24} \left(1^7 + 1.01^7 + \dots + 1.99^7 \right) (0.01)^3 \approx 0.0093731$ ,

which, as expected, is close to the actual error of $102.3 - 102.2904379 \approx 0.00956211$ .

Let $y_i = x_i^{k-2}$ , so that the error becomes

$E \approx \displaystyle \frac{k(k-1)}{24} \left(y_0 + y_1 + \dots + y_{n-1} \right) h^3 + O(h^4) = \displaystyle \frac{k(k-1)}{24} \overline{y} n h^3$ ,

where $\overline{y} = (y_0 + y_1 + \dots + y_{n-1})/n$ is the average of the $y_i$ . Clearly, this average is somewhere between the smallest and the largest of the $y_i$ . Since $y = x^{k-2}$ is a continuous function, that means that there must be some value of $x_*$ between $x_0$ and $x_{k-1}$ — and therefore between $a$ and $b$ — so that $x_*^{k-2} = \overline{y}$ by the Intermediate Value Theorem. We conclude that the error can be written as

$E \approx \displaystyle \frac{k(k-1)}{24} x_*^{k-2} nh^3$ ,

Finally, since $h$ is the length of one subinterval, we see that $nh = b-a$ is the total length of the interval $[a,b]$ . Therefore,

$E \approx \displaystyle \frac{k(k-1)}{24} x_*^{k-2} (b-a)h^2 \equiv ch^2$ ,

where the constant $c$ is determined by $a$ , $b$ , and $k$ . In other words, for the special case $f(x) = x^k$ , we have established that the error from the Midpoint Rule is approximately quadratic in $h$ — without resorting to the generalized mean-value theorem and confirming the numerical observations we made earlier.

Thoughts on Numerical Integration (Part 13): Left endpoint rule and global rate of convergence

Numerical integration is a standard topic in first-semester calculus. From time to time, I have received questions from students on various aspects of this topic, including:

Why is numerical integration necessary in the first place?
Where do these formulas come from (especially Simpson’s Rule)?
How can I do all of these formulas quickly?
Is there a reason why the Midpoint Rule is better than the Trapezoid Rule?
Is there a reason why both the Midpoint Rule and the Trapezoid Rule converge quadratically?
Is there a reason why Simpson’s Rule converges like the fourth power of the number of subintervals?

In the previous post in this series, we found that the local error of the left endpoint approximation to $\displaystyle \int_{x_i}^{x_i+h} x^k \, dx$ was equal to

$\displaystyle \frac{k}{2} x_i^{k-1} h^2 + O(h^3)$ .

The total error when approximating $\displaystyle \int_a^b x^k \, dx = \int_{x_0}^{x_n} x^k , dx$ will be the sum of the errors for the integrals over $[x_0,x_1]$ , $[x_1,x_2]$ , through $[x_{n-1},x_n]$ . Therefore, the total error will be

$E \approx \displaystyle \frac{k}{2} \left(x_0^{k-1} + x_1^{k-1} + \dots + x_{n-1}^{k-1} \right) h^2$ .

$\displaystyle \frac{9}{2} \left(1^8 + 1.01^8 + \dots + 1.99^8 \right) (0.01)^2 \approx 2.49801$ ,

which, as expected, is close to the actual error of $102.3 - 99.76412456 \approx 2.53588$ .

We now perform a more detailed analysis of the global error. Let $y_i = x_i^{k-1}$ , so that the error becomes

$E \approx \displaystyle \frac{k}{2} \left(y_0 + y_1 + \dots + y_{n-1} \right) h^2 + O(h^3) = \displaystyle \frac{k}{2} \overline{y} n h^2$ ,

where $\overline{y} = (y_0 + y_1 + \dots + y_{n-1})/n$ is the average of the $y_i$ . Clearly, this average is somewhere between the smallest and the largest of the $y_i$ . Since $y = x^{k-1}$ is a continuous function, that means that there must be some value of $x_*$ between $x_0$ and $x_{n-1}$ — and therefore between $a$ and $b$ — so that $x_*^{k-1} = \overline{y}$ by the Intermediate Value Theorem. We conclude that the error can be written as

$E \approx \displaystyle \frac{k}{2} x_*^{k-1} nh^2$ ,

Finally, since $h$ is the length of one subinterval, we see that $nh = b-a$ is the total length of the interval $[a,b]$ . Therefore,

$E \approx \displaystyle \frac{k}{2} x_*^{k-1} (b-a)h \equiv ch$ ,

where the constant $c$ is determined by $a$ , $b$ , and $k$ . In other words, for the special case $f(x) = x^k$ , we have established that the error from the left-endpoint rule is approximately linear in $h$ — without resorting to the generalized mean-value theorem.

Thoughts on Numerical Integration (Part 12): Left endpoint rule and local rate of convergence

Numerical integration is a standard topic in first-semester calculus. From time to time, I have received questions from students on various aspects of this topic, including:

Why is numerical integration necessary in the first place?
Where do these formulas come from (especially Simpson’s Rule)?
How can I do all of these formulas quickly?
Is there a reason why the Midpoint Rule is better than the Trapezoid Rule?
Is there a reason why both the Midpoint Rule and the Trapezoid Rule converge quadratically?
Is there a reason why Simpson’s Rule converges like the fourth power of the number of subintervals?

In this post, we will perform an error analysis for the left-endpoint rule

$\int_a^b f(x) \, dx \approx h \left[f(x_0) + f(x_1) + \dots + f(x_{n-1}) \right] \equiv L_n$

where $n$ is the number of subintervals and $h = (b-a)/n$ is the width of each subinterval, so that $x_k = x_0 + kh$ .

As noted above, a true exploration of error analysis requires the generalized mean-value theorem, which perhaps a bit much for a talented high school student learning about this technique for the first time. That said, the ideas behind the proof are accessible to high school students, using only ideas from the secondary curriculum, if we restrict our attention to the special case $f(x) = x^k$ , where $k \ge 5$ is a positive integer.

For this special case, the true area under the curve $f(x) = x^k$ on the subinterval $[x_i, x_i +h]$ will be

$\displaystyle \int_{x_i}^{x_i+h} x^k \, dx = \frac{1}{k+1} \left[ (x_i+h)^{k+1} - x_i^{k+1} \right]$

$= \displaystyle \frac{1}{k+1} \left[x_i^{k+1} + {k+1 \choose 1} x_i^k h + {k+1 \choose 2} x_i^{k-1} h^2 + O(h^3) - x_i^{k+1} \right]$

$= \displaystyle \frac{1}{k+1} \left[ (k+1) x_i^k h + \frac{(k+1)k}{2} x_i^{k-1} h^2 + O(h^3) \right]$

$= x_i^k h + \displaystyle \frac{k}{2} x_i^{k-1} h^2 + O(h^3)$

In the above, the shorthand $O(h^3)$ can be formally defined, but here we’ll just take it to mean “terms that have a factor of $h^3$ or higher that we’re too lazy to write out.” Since $h$ is supposed to be a small number, these terms will be much smaller in magnitude that the terms that have $h$ or $h^2$ and thus can be safely ignored.

Using only the left-endpoint of the subinterval, the left-endpoint approximation of $\displaystyle \int_{x_i}^{x_i+h} x^k \, dx$ is $x_i^k h$ . Therefore, the error in this approximation will be equal to

$\displaystyle \frac{k}{2} x_i^{k-1} h^2 + O(h^3)$ .

In the next post of this series, we’ll show that the global error when integrating between $a$ and $b$ — as opposed to between $x_i$ and $x_i + h$ — is approximately linear in $h$ .

Thoughts on Numerical Integration (Part 11): Simpson’s Rule and exploration of error analysis

Numerical integration is a standard topic in first-semester calculus. From time to time, I have received questions from students on various aspects of this topic, including:

Why is numerical integration necessary in the first place?
Where do these formulas come from (especially Simpson’s Rule)?
How can I do all of these formulas quickly?
Is there a reason why the Midpoint Rule is better than the Trapezoid Rule?
Is there a reason why both the Midpoint Rule and the Trapezoid Rule converge quadratically?
Is there a reason why Simpson’s Rule converges like the fourth power of the number of subintervals?

In the previous post in this series, I discussed three different ways of numerically approximating the definite integral $\displaystyle \int_a^b f(x) \, dx$ , the area under a curve $f(x)$ between $x=a$ and $x=b$ .

In this series, we’ll choose equal-sized subintervals of the interval $[a,b]$ . If $h = (b-a)/n$ is the width of each subinterval so that $x_k = x_0 + kh$ , then the integral may be approximated as

$\int_a^b f(x) \, dx \approx h \left[f(x_0) + f(x_1) + \dots + f(x_{n-1}) \right] \equiv L_n$

using left endpoints,

$\int_a^b f(x) \, dx \approx h \left[f(x_1) + f(x_2) + \dots + f(x_n) \right] \equiv R_n$

using right endpoints, and

$\int_a^b f(x) \, dx \approx h \left[f(c_1) + f(c_2) + \dots + f(c_n) \right] \equiv M_n$

using the midpoints of the subintervals. We have also derived the Trapezoid Rule

$\int_a^b f(x) \, dx \approx \displaystyle \frac{h}{2} [f(x_0) + 2f(x_1) + \dots + 2f(x_{n-1}) + f(x_n)] \equiv T_n$

and Simpson’s Rule (if $n$ is even)

$\int_a^b f(x) \, dx \approx \displaystyle \frac{h}{3} \left[y_0 + 4 y_1 + 2 y_2 + 4 y_3 + \dots + 2y_{n-2} + 4 y_{n-1} + y_{n} \right] \equiv S_n$ .

In the previous post in this series, we saw that both the left-endpoint and right-endpoint rules have a linear rate of convergence: if twice as many subintervals are taken, then the error appears to go down by a factor of 2. If ten times as many subintervals are used, then the error should go down by a factor of 10. However, both the Midpoint Rule and the Trapezoid Rule have a quadratic rate of convergence: if twice as many subintervals are taken, then the error appears to go down by a factor of 4. If ten times as many subintervals are used, then the error should go down by a factor of 100. Moreover, it appears that the error from the Midpoint Rule is about half that of the Trapezoid Rule if the same number of subintervals are used.

Let’s now explore the results of Simpson’s Rule applied to $\displaystyle \int_1^2 x^9 \, dx = 102.3$ using different numbers of subintervals. The results are summarized in the table below.

The first immediate observation is that these approximations are far better than even the Midpoint and Trapezoid Rules! Indeed, we see that $S_{20} \approx 102.301$ , using only 20 subintervals, is a better approximation than (from previous posts) either $M_{100} = 102.290$ or $T_{100} \approx 102.319$ using 100 subintervals! The lesson to learn again: Simpson’s Rule is a bit more difficult to compute than the Trapezoid Rule or the Midpoint Rule because of the different weights. Nevertheless, choosing a good algorithm is often far better than simply doing lots of computations.

There’s a second observation: the rate of convergence appears to be much, much faster. Indeed, the data points appear to fit a quartic polynomial very well. Said another way, if twice as many subintervals are taken, then the error appears to go down by a factor of 16. We can actually see this in the figure, looking at the lines with 10, 20, 40, and 80 subintervals.

Error with 10 subintervals = $|102.3174904 - 102.3| = 0.0174904$ .

Error with 20 subintervals = $|102.3011002 - 102.3| = 0.0011002$ .

Error with 40 subintervals = $|102.3000689 - 102.3| = 0.0000689$ .

Error with 80 subintervals = $|102.3000043 - 102.3| = 0.0000043$ .

In all cases, the error on the next line is about the error on the previous line divided by 16.

This illustrates quartic convergence, as opposed to the mere linear convergence of the left- and right-endpoint rules or the quadratic convergence of the Midpoint and Trapezoid Rules.

Thoughts on Numerical Integration (Part 10): Trapezoid Rule and exploration of error analysis

Numerical integration is a standard topic in first-semester calculus. From time to time, I have received questions from students on various aspects of this topic, including:

Why is numerical integration necessary in the first place?
Where do these formulas come from (especially Simpson’s Rule)?
How can I do all of these formulas quickly?
Is there a reason why the Midpoint Rule is better than the Trapezoid Rule?
Is there a reason why both the Midpoint Rule and the Trapezoid Rule converge quadratically?
Is there a reason why Simpson’s Rule converges like the fourth power of the number of subintervals?

In this series, we’ll choose equal-sized subintervals of the interval $[a,b]$ . If $h = (b-a)/n$ is the width of each subinterval so that $x_k = x_0 + kh$ , then the integral may be approximated as

$\int_a^b f(x) \, dx \approx h \left[f(x_0) + f(x_1) + \dots + f(x_{n-1}) \right] \equiv L_n$

using left endpoints,

$\int_a^b f(x) \, dx \approx h \left[f(x_1) + f(x_2) + \dots + f(x_n) \right] \equiv R_n$

using right endpoints, and

$\int_a^b f(x) \, dx \approx h \left[f(c_1) + f(c_2) + \dots + f(c_n) \right] \equiv M_n$

using the midpoints of the subintervals. We have also derived the Trapezoid Rule

$\int_a^b f(x) \, dx \approx \displaystyle \frac{h}{2} [f(x_0) + 2f(x_1) + \dots + 2f(x_{n-1}) + f(x_n)] \equiv T_n$

and Simpson’s Rule (if $n$ is even)

$\int_a^b f(x) \, dx \approx \displaystyle \frac{h}{3} \left[y_0 + 4 y_1 + 2 y_2 + 4 y_3 + \dots + 2y_{n-2} + 4 y_{n-1} + y_{n} \right] \equiv S_n$ .

Let’s now explore the results of the Trapezoid Rule applied to $\displaystyle \int_1^2 x^9 \, dx = 102.3$ using different numbers of subintervals. The results are summarized in the table below.

Once again, the data points fit a quadratic polynomial well, indicating quadratic convergence.

More subtly, it appears that the Trapezoid Rule isn’t quite as good as the Midpoint Rule. Here are the results from the Midpoint Rule (which also appeared in the previous post in this series):

For $n = 100$ subintervals, the error of the Trapezoid Rule is $|102.3191-102.3| = 0.0191$ , which the error from the Midpoint Rule is $|102.2904-102.3| = 0.0096$ . In other words, while both of these methods are superior to the left- and right-endpoint rules, it appears that the error from the Midpoint Rule is about half of the error from the Trapezoid Rule. The Midpoint Rule appears to be better.

To me, this is far from an obvious conclusion. Geometrically, it’s far from clear that the rectangles from the Midpoint Rule…

… provide a better approximation than using trapezoids …

… yet it appears that’s exactly what happened. This can be rigorously proven, as we’ll explore later in this series.

Thoughts on Numerical Integration (Part 9): Midpoint rule and exploration of error analysis

Numerical integration is a standard topic in first-semester calculus. From time to time, I have received questions from students on various aspects of this topic, including:

Why is numerical integration necessary in the first place?
Where do these formulas come from (especially Simpson’s Rule)?
How can I do all of these formulas quickly?
Is there a reason why the Midpoint Rule is better than the Trapezoid Rule?
Is there a reason why both the Midpoint Rule and the Trapezoid Rule converge quadratically?
Is there a reason why Simpson’s Rule converges like the fourth power of the number of subintervals?

In this series, we’ll choose equal-sized subintervals of the interval $[a,b]$ . If $h = (b-a)/n$ is the width of each subinterval so that $x_k = x_0 + kh$ , then the integral may be approximated as

$\int_a^b f(x) \, dx \approx h \left[f(x_0) + f(x_1) + \dots + f(x_{n-1}) \right] \equiv L_n$

using left endpoints,

$\int_a^b f(x) \, dx \approx h \left[f(x_1) + f(x_2) + \dots + f(x_n) \right] \equiv R_n$

using right endpoints, and

$\int_a^b f(x) \, dx \approx h \left[f(c_1) + f(c_2) + \dots + f(c_n) \right] \equiv M_n$

using the midpoints of the subintervals. We have also derived the Trapezoid Rule

$\int_a^b f(x) \, dx \approx \displaystyle \frac{h}{2} [f(x_0) + 2f(x_1) + \dots + 2f(x_{n-1}) + f(x_n)] \equiv T_n$

and Simpson’s Rule (if $n$ is even)

$\int_a^b f(x) \, dx \approx \displaystyle \frac{h}{3} \left[y_0 + 4 y_1 + 2 y_2 + 4 y_3 + \dots + 2y_{n-2} + 4 y_{n-1} + y_{n} \right] \equiv S_n$ .

In the previous post in this series, we saw that both the left-endpoint and right-endpoint rules have a linear rate of convergence: if twice as many subintervals are taken, then the error appears to go down by a factor of 2. If ten times as many subintervals are used, then the error should go down by a factor of 10. Let’s now explore the results of the midpoint rule applied to $\displaystyle \int_1^2 x^9 \, dx = 102.3$ using different numbers of subintervals. The results are summarized in the table below. The first immediate observation is that these approximations are far better than the left- and right-endpoint rule approximations! Indeed, we see that $M_{10} \approx 101.3$ , using only ten subintervals, is a far better approximation than (from the previous post) either $L_{100} = 99.8$ or $R_{100} \approx 104.9$ using 100 subintervals! The lesson to learn: choosing a good algorithm is often far better than simply doing lots of computations.

There’s a second observation: the rate of convergence appears to be much, much faster. Indeed, the data points appear to fit a parabola very well instead of a straight line. Said another way, if twice as many subintervals are taken, then the error appears to go down by a factor of 4. If ten times as many subintervals are used, then the error should go down by a factor of 100. This illustrates quadratic convergence, as opposed to the mere linear convergence of the left- and right-endpoint rules.

Thoughts on Numerical Integration (Part 8): Left and right endpoint rules and exploration of error analysis

Numerical integration is a standard topic in first-semester calculus. From time to time, I have received questions from students on various aspects of this topic, including:

Why is numerical integration necessary in the first place?
Where do these formulas come from (especially Simpson’s Rule)?
How can I do all of these formulas quickly?
Is there a reason why the Midpoint Rule is better than the Trapezoid Rule?
Is there a reason why both the Midpoint Rule and the Trapezoid Rule converge quadratically?
Is there a reason why Simpson’s Rule converges like the fourth power of the number of subintervals?

In this series, we’ll choose equal-sized subintervals of the interval $[a,b]$ . If $h = (b-a)/n$ is the width of each subinterval so that $x_k = x_0 + kh$ , then the integral may be approximated as

$\int_a^b f(x) \, dx \approx h \left[f(x_0) + f(x_1) + \dots + f(x_{n-1}) \right] \equiv L_n$

using left endpoints,

$\int_a^b f(x) \, dx \approx h \left[f(x_1) + f(x_2) + \dots + f(x_n) \right] \equiv R_n$

using right endpoints, and

$\int_a^b f(x) \, dx \approx h \left[f(c_1) + f(c_2) + \dots + f(c_n) \right] \equiv M_n$

using the midpoints of the subintervals. We have also derived the Trapezoid Rule

$\int_a^b f(x) \, dx \approx \displaystyle \frac{h}{2} [f(x_0) + 2f(x_1) + \dots + 2f(x_{n-1}) + f(x_n)] \equiv T_n$

and Simpson’s Rule (if $n$ is even)

$\int_a^b f(x) \, dx \approx \displaystyle \frac{h}{3} \left[y_0 + 4 y_1 + 2 y_2 + 4 y_3 + \dots + 2y_{n-2} + 4 y_{n-1} + y_{n} \right] \equiv S_n$ .

All of the above approximations to $\displaystyle \int_a^b f(x) \, dx$ are precisely that — approximations. That begs the obvious question: how can we get better approximations. One obvious answer is taking more subintervals. The figures below show the left-endpoint approximations using $n = 4$ and $n = 40$ subintervals. Geometrically, it’s clear that the orange rectangles in the second picture do a better job of approximating the area under the curve. Unfortunately, simply taking more subintervals has its limitations. Using a spreadsheet as in the previous post in this series, one can implement 100 or even 1000 subintervals without much difficult. However, as demonstrated in the video below, implementing any of these methods with 10,000 subintervals is pretty time-consuming. (Tl/dw: It can take literally a couple of minutes.)

Instead of relying on sheer computational firepower, let’s instead investigate how good these numerical methods actually are. To begin, let’s explore the left-endpoint rule applied to $\displaystyle \int_1^2 x^9 \, dx = 102.3$ using different numbers of subintervals. The results are summarized in the table below. As $n$ increases, the left endpoint approximations $L_n$ are indeed getting closer and closer to the actual value of $\displaystyle \int_1^2 x^9 \, dx = 102.3$ . Interestingly, when the width $h$ is plotted with these approximations, the data points fall almost exactly on a straight line. The same phenomenon occurs when using right endpoints: So, it appears that the errors in both the left- and right-endpoint rules are a linear function of the size of the subintervals. Said another way, if twice as many subintervals are taken, then the error appears to go down by a factor of 2. If ten times as many subintervals are used, then the error should go down by a factor of 10. As we’ll see in the next few posts, the errors for the Midpoint Rule, the Trapezoid Rule, and especially Simpson’s Rule are much better than the errors from these two methods.

Thoughts on Numerical Integration (Part 7): Implementation with Excel

Numerical integration is a standard topic in first-semester calculus. From time to time, I have received questions from students on various aspects of this topic, including:

Why is numerical integration necessary in the first place?
Where do these formulas come from (especially Simpson’s Rule)?
How can I do all of these formulas quickly?
Is there a reason why the Midpoint Rule is better than the Trapezoid Rule?
Is there a reason why both the Midpoint Rule and the Trapezoid Rule converge quadratically?
Is there a reason why Simpson’s Rule converges like the fourth power of the number of subintervals?

In this series, we’ll choose equal-sized subintervals of the interval $[a,b]$ . If $h = (b-a)/n$ is the width of each subinterval so that $x_k = x_0 + kh$ , then the integral may be approximated as

$\int_a^b f(x) \, dx \approx h \left[f(x_0) + f(x_1) + \dots + f(x_{n-1}) \right] \equiv L_n$

using left endpoints,

$\int_a^b f(x) \, dx \approx h \left[f(x_1) + f(x_2) + \dots + f(x_n) \right] \equiv R_n$

using right endpoints, and

$\int_a^b f(x) \, dx \approx h \left[f(c_1) + f(c_2) + \dots + f(c_n) \right] \equiv M_n$

using the midpoints of the subintervals. We have also derived the Trapezoid Rule

$\int_a^b f(x) \, dx \approx \displaystyle \frac{h}{2} [f(x_0) + 2f(x_1) + \dots + 2f(x_{n-1}) + f(x_n)] \equiv T_n$

and Simpson’s Rule (if $n$ is even)

$\int_a^b f(x) \, dx \approx \displaystyle \frac{h}{3} \left[y_0 + 4 y_1 + 2 y_2 + 4 y_3 + \dots + 2y_{n-2} + 4 y_{n-1} + y_{n} \right] \equiv S_n$ .

Computing any of the above formulas on a hand-held calculator can tax the patience of even the most error-conscious student. Indeed, I prefer that my students, when first learning these concepts, use a spreadsheet instead of a calculator or even a computer program, as I think that the visual layout of the spreadsheet aids in understanding how the formula works. In what follows, I implement the above formulas for the integral $\displaystyle \int_1^2 x^9 \, dx$ using $n=10$ subintervals, so that $h = (2-1)/10 = 0.1$ . To implement the left-endpoint rule, I enter the labels “x” and “x^9” in cells A1 and B1 of a spreadsheet. I then enter 1 (the left endpoint) in cell A2. In cell A3, I enter “=A2+0.1”, instructing the spreadsheet to add 0.1 to the value in cell A2. Then, instead of typing all of the other values of $x_k$ , I use the fill-down feature to repeat this pattern for cells A3 through A11. In cell B2, I enter “=A1^9”, applying the function $f(x) = x^9$ to the $x-$ coordinate in cell A2. Again, I use the fill-down feature to repeat this pattern for cells B3-B11. The fill-down feature saves a lot of time! Finally, in cell B13, I enter “=0.1*SUM(B2:B11)”, adding the values in cells B2 through B11 and multiplying the sum by $h$ . The result, $78.6581$ , is the approximation using the left-endpoint rule with 10 subintervals. Once this is done, the right-endpoint rule can be obtained almost for free. The only change is to change the value of cell A2 from 1 to 1.1. Everything else should automatically update. The midpoint rule is also obtained quickly by changing the value of cell A2 from 1 to 1.05, the midpoint of the first subinterval $[1,1.01]$ . Implementing the Trapezoid Rule requires a little more work. We reset the value of A2 back to 1, the value of the left-endpoint. We also fill down the pattern one extra row (in this case, row 12). To implement the Trapezoid Rule, we have to multiply all function values (except for those at the endpoints) by 2. To implement this, I introduce column C. These weights can be typed by hand, but again the fill-down feature can speed things up. Then, in column D, I multiply the values in columns B and C. For example, the result in cell D2 is obtained by typing “=B2*C2”. Once again, the fill-down feature is used for all rows. Finally, the approximation itself is obtained by typing “=0.1/2*SUM(D2:D12)” in cell D13. After implementing the Trapezoid Rule, Simpson’s Rule is not much more effort. The biggest change is the alternating weights, so that the endpoints have weight 1 while the others oscillate between 4 and 2, ending on 4 on the second-to-last value of $x$ . Again, these could be typed by hand, but it’s easiest to enter 4 in cell C3, 2 in cell C4, and then “=C3” in cell C5. The fill-down feature can take care of the rest of the weights. The Simpson’s Rule approximation is obtained by typing “=0.1/3*SUM(D2:D12)” in cell D13, with a new denominator of 3.

Thoughts on Numerical Integration (Part 6): Connection between Simpson’s Rule, Trapezoid Rule, and Midpoint Rule

Numerical integration is a standard topic in first-semester calculus. From time to time, I have received questions from students on various aspects of this topic, including:

Why is numerical integration necessary in the first place?
Where do these formulas come from (especially Simpson’s Rule)?
How can I do all of these formulas quickly?
Is there a reason why the Midpoint Rule is better than the Trapezoid Rule?
Is there a reason why both the Midpoint Rule and the Trapezoid Rule converge quadratically?
Is there a reason why Simpson’s Rule converges like the fourth power of the number of subintervals?

In this series, we’ll choose equal-sized subintervals of the interval $[a,b]$ . If $h = (b-a)/n$ is the width of each subinterval so that $x_k = x_0 + kh$ , then the integral may be approximated as

$\int_a^b f(x) \, dx \approx h \left[f(x_0) + f(x_1) + \dots + f(x_{n-1}) \right] \equiv L_n$

using left endpoints,

$\int_a^b f(x) \, dx \approx h \left[f(x_1) + f(x_2) + \dots + f(x_n) \right] \equiv R_n$

using right endpoints, and

$\int_a^b f(x) \, dx \approx h \left[f(c_1) + f(c_2) + \dots + f(c_n) \right] \equiv M_n$

using the midpoints of the subintervals. We have also derived the Trapezoid Rule

$\int_a^b f(x) \, dx \approx \displaystyle \frac{h}{2} [f(x_0) + 2f(x_1) + \dots + 2f(x_{n-1}) + f(x_n)] \equiv T_n$

and Simpson’s Rule (if $n$ is even)

$\int_a^b f(x) \, dx \approx \displaystyle \frac{h}{3} \left[y_0 + 4 y_1 + 2 y_2 + 4 y_3 + \dots + 2y_{n-2} + 4 y_{n-1} + y_{n} \right] \equiv S_n$ .

There is a somewhat surprising connection between the last three formulas. Let’s divide the interval $[a,b]$ into $2n$ subintervals with $h = (b-a)/(2n)$ and $x_0 = a$ , $x_1 = x_0 + h$ , $x_2 = x_0 + 2h$ , and so on. Then Simpson’s Rule becomes

$S_{2n} = \displaystyle \frac{h}{3} \left[y_0 + 4 y_1 + 2 y_2 + 4 y_3 + \dots + 2y_{2n-2} + 4 y_{2n-1} + y_{2n} \right]$ .

Next, let’s divide the interval $[a,b]$ into $n$ subintervals, but let’s not redefine the values of $h$ and the $x_k$ . Instead, the width of each subinterval will be $(b-a)/n$ , which is equal to $2h$ . (In other words, since there are half as many subintervals, each one is twice as long.) Also, the endpoints of these subintervals will be $x_0 = a$ , $x_2 = x_0 + 2h$ , $x_4 = x_0 + 4h$ , and so on. So, keeping the same labeling convention as with Simpson’s Rule, the Trapezoid Rule becomes

$T_n = \displaystyle \frac{2h}{2} [f(x_0) + 2f(x_2) + 2f(x_4) + \dots + 2f(x_{2n-2}) + f(x_{2n})]$

$= h [f(x_0) + 2f(x_2) + 2f(x_4) + \dots + 2f(x_{2n-2}) + f(x_{2n})]$ .

(Again, the width of the subintervals in this case is $2h$ , where $h = (b-a)/2n$ .) Furthermore, the midpoint of subinterval $[x_0, x_2]$ will be $x_1$ , the midpoint of subinterval $[x_2,x_4]$ will be $x_3$ , and so on. Therefore, keeping the same labeling convention, the Midpoint Rule becomes

$M_n = \displaystyle 2h [f(x_1) + f(x_3) + f(x_5) + \dots + f(x_{2n-1}) ]$ .

It turns out that $\displaystyle \frac{2}{3} M_n + \frac{1}{3} T_n$ , a certain weighted average of $T_n$ and $M_n$ , is equal to

$\displaystyle \frac{4h}{3} [f(x_1) + f(x_3) + \dots + f(x_{2n-1}) ] + \frac{h}{3} [f(x_0) + 2f(x_2) + \dots + 2f(x_{2n-2}) + f(x_{2n})]$

$= \displaystyle \frac{h}{3} [f(x_0) + 4 f(x_1) + 2f(x_2) + \dots + 2f(x_{2n-2}) + 4 f(x_{2n-1} + f(x_{2n})]$

$= S_{2n}$ .

So, if the Midpoint Rule and the Trapezoid Rule have already been computed for $n$ subintervals, then Simpson’s Rule for $2n$ subintervals can be computed at almost no additional effort.

Thoughts on Numerical Integration (Part 4): Derivation of Trapezoid Rule

Numerical integration is a standard topic in first-semester calculus. From time to time, I have received questions from students on various aspects of this topic, including:

Why is numerical integration necessary in the first place?
Where do these formulas come from (especially Simpson’s Rule)?
How can I do all of these formulas quickly?
Is there a reason why the Midpoint Rule is better than the Trapezoid Rule?
Is there a reason why both the Midpoint Rule and the Trapezoid Rule converge quadratically?
Is there a reason why Simpson’s Rule converges like the fourth power of the number of subintervals?

In this series, we’ll choose equal-sized subintervals of the interval $[a,b]$ . If $h = (b-a)/n$ is the width of each subinterval so that $x_k = x_0 + kh$ , then the integral may be approximated as

$\int_a^b f(x) \, dx \approx h \left[f(x_0) + f(x_1) + \dots + f(x_{n-1}) \right] \equiv L_n$

using left endpoints,

$\int_a^b f(x) \, dx \approx h \left[f(x_1) + f(x_2) + \dots + f(x_n) \right] \equiv R_n$

using right endpoints, and

$\int_a^b f(x) \, dx \approx h \left[f(c_1) + f(c_2) + \dots + f(c_n) \right] \equiv M_n$

using the midpoints of the subintervals.

All three of these approximations were obtained by approximating the above shaded region by rectangles. However, perhaps it might be better to use some other shape besides rectangles. In the Trapezoidal Rule, we approximate the area by using (surprise!) trapezoids, as in the figure below.

The first trapezoid has height $h$ and bases $f(x_0)$ and $f(x_1)$ , and so the area of the first trapezoid is $\frac{1}{2} h[ f(x_0) + f(x_1) ]$ . The other areas are found similarly. Adding these together, we get the approximation

$T_n = \displaystyle \frac{h}{2}[f(x_0) + f(x_1)] + \frac{h}{2} [f(x_1) + f(x_2)] + \dots +$

$+ \displaystyle \frac{h}{2} [f(x_{n-2})+f(x_{n-1})] + \frac{h}{2} [f(x_{n-1})+f(x_n)]$

$= \displaystyle \frac{h}{2} [f(x_0) + 2f(x_1) + 2f(x_2) + \dots + 2f(x_{n-2}) + 2f(x_{n-1}) + f(x_n)].$

Interestingly, $T_n$ is the average of the two endpoint approximations $L_n$ and $R_n$ :

$\displaystyle \frac{L_n+R_n}{2} = \frac{L_n}{2} + \frac{R_n}{2}$

$= \displaystyle \frac{h}{2} \left[f(x_0) + f(x_1) + f(x_2) + \dots + f(x_{n-1}) \right]$

$+\displaystyle \frac{h}{2} \left[f(x_1) + f(x_2) + \dots + f(x_{n-1}) + f(x_{n}) \right]$

$= \displaystyle \frac{h}{2} \left[f(x_0) + 2f(x_1) + \dots + 2f(x_{n-1}) + f(x_n) \right]$

$= T_n$ .

Of course, as a matter of computation, it’s a lot quicker to directly compute $T_n$ instead of computing $L_n$ and $R_n$ separately and then averaging.