This section concentrates on definite integrals. The inputs are y(x) and two endpoints a and b. The output is the integral I. Our goal is to find that number ∫ab y(x) dx = I, accurately and in a short time. Normally this goal is achievable—as soon as we have a good method for computing integrals.
Our two approaches so far have weaknesses. The search for an antiderivative succeeds in important cases, and Chapter 7 extends that range—but generally ƒ(x) is not available. The other approach (by rectangles) is in the right direction but too crude. The height is set by y(x) at the right and left end of each small interval. The right and left rectangle rules add the areas (Δx times y):
Rn = (Δx)(y1 + y2 + … + yn) and Ln = (Δx)(y0 + y1 + … + yn-1).
The value of y(x) at the end of interval j is yj. The extreme left value yo = y(a) enters Ln. With n equal intervals of length Δx = (b−a)/n, the extreme right value is yn = y(b). It enters Rn. Otherwise the sums are the same—simple to compute, easy to visualize, but very inaccurate.
This section goes from slow methods (rectangles) to better methods (trapezoidal and midpoint) to good methods (Simpson and Gauss). Each improvement cuts down the error. You could discover the formulas without the book, by integrating x and x² and x⁴. The rule Rn would come out on one side of the answer, and Ln would be on the other side. You would figure out what to do next, to come closer to the exact integral. The book can emphasize one key point:
The quality of a formula depends on how many integrals ∫ 1 dx, ∫ x dx, ∫ x² dx, …, it computes exactly. If ∫ xp dx is the first to be wrong, the order of accuracy is p.
By testing the integrals of 1, x, x², …, we decide how accurate the formulas are.
Figure 5.18 shows the rectangle rules Rn and Ln. They are already wrong when y = x. These methods are first-order: p = 1. The errors involve the first power of Δx—where we would much prefer a higher power. A larger p in (Δx)p means a smaller error.
When the graph of y(x) is a straight line, the integral I is known. The error triangles E and e have base Δx. Their heights are the differences yj+1 − yj. The areas are ½(base)(height), and the only difference is a minus sign. (L is too low, so the error L − I is negative.) The total error in Rn is the sum of the E's:
Rn − I = ½Δx(y1 - Y0) + … + ½Δx(yn − yn-1)= ½Δx(yn − y0). &enmsp; (1)
All y's between y0 and yn cancel. Similarly for the sum of the e's:
Ln − I = −½Δx(yn − y0) = −½Δx[y(b) − y(a)]. (2)
The greater the slope of y(x), the greater the error—since rectangles have zero slope.
Formulas (1) and (2) are nice—but those errors are large. To integrate y = x from a = 0 to b = 1, the error is &frac1;Δx(1 − 0). It takes 500,000 rectangles to reduce this error to 1/1,000,000. This accuracy is reasonable, but that many rectangles is unacceptable.
The beauty of the error formulas is that they are "asymptotically correct" for all functions. When the graph is curved, the errors don't fit exactly into triangles. But the ratio of predicted error to actual error approaches 1. As Δx → 0, the graph is almost straight in each interval—this is linear approximation.
The error prediction ½Δx[y(b) − y(a)] is so simple that we test it on y(x) = √x:
| I = | ∫ | 1 | √x dx = 2/3 |
| 0 |
| n = | 1 | 10 | 100 | 1000 |
| error Rn − I = | .33 | .044 | .0048 | .00049 |
| error Ln − I = | −.67 | −.056 | −.0052 | −.00051 |
The error decreases along each row. So does Δx = .1, .01, .001, .0001. Multiplying n by 10 divides Δx by 10. The error is also divided by 10 (almost). The error is nearly proportional to Δx—typical of first-order methods.
The predicted error is ½Δx, since here y(1) = 1 and y(0) = 0. The computed errors in the table come closer and closer to ½Δx = .5, .05, .005, .0005. The prediction is the "leading term" in the actual error.
The table also shows a curious fact. Subtracting the last row from the row above gives exact numbers 1, .1, .01, and .001. This is (Rn − I) − (Ln − I), which is Rn − Ln. It comes from an extra rectangle at the right, included in Rn but not Ln. Its height is 1 and its area is 1, .1, .01, .001.
The errors in Rn and Ln almost cancel. The average Tn = ½(Rn + Ln) has less error—it is the "trapezoidal rule." First we give the rectangle prediction two final tests:
| n=1 | n=10 | n=100 | n=1000 | ||
| ∫ (x² − x) dx: | errors | 1.7⋅10-1 | 1.7⋅10-3 | 1.7⋅10-5 | 1.7⋅10-7 |
| ∫ dx/(10 + cos 2πx): | errors | −1⋅10-3 | 2⋅10-14 | "0" | "0" |
Those errors are falling faster than Δx. For y = x² − x the prediction explains why: y(0) equals y(1). The leading term, with y(b) minus y(a), is zero. The exact errors are (1/6)(Δx)², dropping from 10-1 to 10-3 to 10-5 to 10-7. In these examples Ln is identical to Rn (and also to Tn), because the end rectangles are the same. We will see these (1/6)(Δx)² errors in the trapezoidal rule.
The last row in the table is more unusual. It shows practically no error. Why do the rectangle rules suddenly achieve such an outstanding success?
The reason is that y(x) = 1/(10 + cos 2πx) is periodic. The leading term in the error is zero, because y(0) = y(1). Also the next term will be zero, because y'(0) = y'(1). Every power of Δx is multiplied by zero, when we integrate over a complete period. So the errors go to zero exponentially fast.
Personal note I tried to integrate 1/(10 + cos 2πx) by hand and failed. Then I was embarrassed to discover the answer in my book on applied mathematics. The method was a special trick using complex numbers, which applies over an exact period. Finally I found the antiderivative (quite complicated) in a handbook of integrals, and verified the area 1/√99.
We move to integration formulas that are exact when y = x. They have second-order accuracy. The Δx error term disappears. The formulas give the correct area under straight lines. The predicted error is a multiple of (Δx)². That multiple is found by testing y = x²—for which the answers are not exact.
The first formula combines Rn and Ln. To cancel as much error as possible, take the average ½(Rn + Ln). This yields the trapezoidal rule, which approximates ∫ y(x) dx by Tn:
Tn = ½Rn + ½Ln = Δ(½y0 + y1 + y2 + … + yn-1 + ½yn). (3)
Another way to find Tn is from the area of the "trapezoid" below y = x in Figure 5.19a.
The base is Δx and the sides have heights yj-1 and yj. Adding those areas gives ½(Ln + Rn) in formula (3)—the coefficients of yj combine into ½ + ½ = 1. Only the first and last intervals are missing a neighbor, so the rule has ½y0 and ½yn. Because trapezoids (unlike rectangles) fit under a sloping line, Tn is exact when y = x.
What is the difference from rectangles? The trapezoidal rule gives weight ½Δx to y0 and yn. The rectangle rule Rn gives full weight Δx to yn (and no weight to y0). Rn − Tn is exactly the leading error ½yn − ½y0. The change to Tn knocks out that error.
Another important formula is exact for y(x) = x. A rectangle has the same area as a trapezoid, if the height of the rectangle is halfway between yj-1 and yj. On a straight line graph that is achieved at the midpoint of the interval. By evaluating y(x) at the halfway points 1/2Δx, 3/2Δx, 5/2Δx, …, we get much better rectangles. This leads to the midpoint rule Mn:
Mn = Δx(y1/2 + y3/2 + … + yn-1/2). (4)
For ∫04 x dx, trapezoids give ½(0) + 1 + 2 + 3 + ½(4) = 8. The midpoint rule gives 1/2 + 3/2 + 5/2 + 7/2 = 8, again correct. The rules become different when y = x², because y1/2 is no longer the average of y0 and y1. Try both second-order rules on x²:
| I = | ∫ | 1 | x² dx |
| 0 |
| n = | 1 | 10 | 100 |
| error Tn − I = | 1/6 | 1/600 | 1/60000 |
| error Mn − I = | −1/12 | −1/1200 | −1/120000 |
The errors fall by 100 when n is multiplied by 10. The midpoint rule is twice as good (−1/12 vs. 1/6). Since all smooth functions are close to parabolas (quadratic approximation in each interval), the leading errors come from Figure 5.19. The trapezoidal error is exactly (1/6)(Δx)² when y(x) is x² (the 12 in the formula divides the 2 in y'):
Tn − I ≈ (1/12)(Δx)²[(y'1 − y'0) + … (y'n − y'n-1)] = (1/12)(Δx)²[y'n − y'0] (5)
Mn − I ≈ −(1/24)(Δx)²[y'n − y'0] = −(1/24)(Δx)²[y'(b) − y'(a)] (6)
For exact error formulas, change y'(b) − y'(a) to (b − a)y''(c). The location of c is unknown (as in the Mean Value Theorem). In practice these formulas are not much used—they involve the pth derivative at an unknown location c. The main point about the error is the factor (Δx)p.
One crucial fact is easy to overlook in our tests. Each value of y(x) can be extremely hard to compute. Every time a formula asks for yj, a computer calls a subroutine. The goal of numerical integration is to get below the error tolerance, while calling for a minimum number of values of y. Second-order rules need about a thousand values for a typical tolerance of 10-6. The next methods are better.
The trapezoidal error is nearly twice the midpoint error (1/6 vs. −1/12). So a good combination will have twice as much of Mn as Tn. That is Simpson's rule:
Sn = (1/3)Tn + (2/3)Mn = (1/6)Δx[y0 + 4y1/2 + 2y1 + 4y3/2 + 2y2 + … + 4yn-1/2 + yn]. (7)
Multiply the midpoint values by 2/3 ≈ 4/6. The endpoint values are multiplied by 2/6, except at the far ends a and b (with heights y0 and yn). This 1-4-2-4-2-4-1 pattern has become famous.
Simpson's rule goes deeper than a combination of T and M. It comes from a parabolic approximation to y(x) in each interval. When a parabola goes through y0, y1/2, y1, the area under it is (1/6)Δx(y0 + 4y1/2 + y1). This is S over the first interval. All our rules are constructed this way: Integrate correctly as many powers 1, x, x², … as possible. Parabolas are better than straight lines, which are better than flat pieces. S beats M, which beats R. Check Simpson's rule on powers of x, with Δx = 1/n:
| n = 1 | n = 10 | n = 100 | |
| error y = x² | 0 | 0 | 0 |
| error y = x³ | 0 | 0 | 0 |
| error y = x⁴ | 8.33⋅10-3 | 8.33⋅10-7 | 8.33⋅10-11 |
Exact answers for x² are no surprise. Sn was selected to get parabolas right. But the zero errors for x³ were not expected. The accuracy has jumped to fourth order, with errors proportional to (Δx)⁴. That explains the popularity of Simpson's rule.
To understand why x³ is integrated exactly, look at the interval [-1, 1]. The odd function x³ has zero integral, and Simpson agrees by symmetry:
| ∫ | 1 | x³ dx = | 1 | x⁴ | 1 | = 0 and | 2 | [(-1)³ + 4(0)³ + 1³] = 0. (8) | |
| -1 | 4 | -1 | 6 |
We need a competitor for Simpson, and Gauss can compete with anybody. He calculated integrals in astronomy, and discovered that two points are enough for a fourth-order method. From -1 to 1 (a single interval) his rule is
| ∫ | 1 | y(x) dx ≈ y(-1/√3) + y(1/√3). (8) |
| -1 |
Those "Gauss points" x = -1/√3 and x = 1/√3 can be found directly. By placing them symmetrically, all odd powers x, x³, … are correctly integrated. The key is in y = x², whose integral is 2/3. The Gauss points -xG and +xG get this integral right:
2/3 = (-xG)² + (xG)², so xG² = 1/3 and xG = ±1/√3.
Figure 5.20c shifts to the interval from 0 to Δx. The Gauss points are (1 ± 1/√3)Δx/2. They are not as convenient as Simpson's (which hand calculators prefer). Gauss is good for thousands of integrations over one interval. Simpson is good when intervals go back to back—then Simpson also uses two y's per interval. For y = x⁴, you see both errors drop by 10-4 in comparing n = 1 to n = 10:
| I = | ∫ | 1 | x⁴ dx | Simpson error | 8.33⋅10-3 | 8.33⋅10-7 |
| 0 | Gauss error | -5.56⋅10-3 | -5.56⋅10-7 |
It is fascinating to know how numerical integration is actually done. The points are not equally spaced! For an integral from 0 to 1, Hewlett-Packard machines might internally replace x by 3u² − 2u³ (the limits on u are also 0 and 1). The machine remembers to change dx. For example,
| ∫ | 1 | dx | becomes | ∫ | 1 | 6(u − u²) du | = | ∫ | 1 | 6(1 − u) du |
| 0 | √x | 0 | √3u² − 2u³ | 0 | √3 − 2u |
Algebraically that looks worse—but the infinite value of 1/√x at x = 0 disappears at u = 0. The differential 6(u − u²) du was chosen to vanish at u = 0 and u = 1. We don't need y(x) at the endpoints—where infinity is most common. In the u variable the integration points are equally spaced—therefore in x they are not.
When a difficult point is inside [a, b], break the interval in two pieces. And chop off integrals that go out to infinity. The integral of ex² should be stopped by x = 10, since the tail is so thin. (It is bad to go too far.) Rapid oscillations are among the toughest—the answer depends on cancellation of highs and lows, and the calculator requires many integration points.
The change from x to u affects periodic functions. I thought equal spacing was good, since 1/(10 + cos 2πx) was integrated above to enormous accuracy. But there is a danger called aliasing. If sin 8πx is sampled with Δx = 1/8, it is always zero. A high frequency 8 is confused with a low frequency 0 (its "alias" which agrees at the sample points). With unequal spacing the problem disappears. Notice how any integration method can be deceived:
Ask for the integral of y = 0 and specify the accuracy. The calculator samples y at x1, …, xk. (With a PAUSE key, the x's may be displayed.) Then integrate Y(x) = (x − x1)² … (x − xk)². That also returns the answer zero (now wrong), because the calculator follows the same steps.
On the HP-28s you enter the function, the endpoints, and the accuracy. The variable x can be named or not (see the margin). The outputs 4.67077 and 4.7E-5 are the requested integral ∫12 ex dx and the estimated error bound. Your input accuracy .00001 guarantees
| 3: 'EXP(X)' 2: {X 1 2} 1: .00001 | relative error in y = | true y − computed y | ≤ .00001. | 3: ((EXP)) 2: {1 2} 1: .00001 |
| computed y |
The machine estimates accuracy based on its experience in sampling y(x). If you guarantee ex within .00000000001, it thinks you want high accuracy and takes longer.
In consulting for HP, William Kahan chose formulas using 1, 3, 7, 15, … sample points. Each new formula uses the samples in the previous formula. The calculator stops when answers are close. The last paragraphs are based on Kahan's work.
TI-81 Program to Test the Integration Methods L, R, T, M, S
| PrgmI:NUM INT :Disp "A=" :Input A :Disp "B=" :Input B :Lbl N :Disp "N=" :Input N :(B-A)/N→D |
:D/2→H :A→X :Y1→L :1→J :0→R :0→M :Lbl 1 :X+H→X :M+Y1→M |
:A+JD→X :R+Y1→R :IS>(J,N) :Goto 1 :(L+R−Y1)→L :RD→R :MD→M :(L+R)/2→T :(2M+T)/3→S |
:Disp "L, R, M, T, S" :Disp L :Disp R :Disp M :Disp T :Disp S :Pause :Goto N |
Place the integrand y(x) in the Y 1 position on the Y = function edit screen. Execute
this program, indicating the interval [A, B] and the number of subintervals N. Rules
L and R and M use N evaluations of y(x). The trapezoidal rule uses N+1 and Simpson's rule uses 2N +1. The program pauses to display the results. Press ENTER to continue by choosing a different N. The program never terminates (only pauses). You break out by pressing ON. Don't forget that IS, Goto, … are on menus.
