The integral

measures an area under the bell curve, and it cannot be evaluated by the fundamental theorem of calculus in the usual way: \(e^{-x^2}\) has no antiderivative that can be written with elementary functions. The number is perfectly definite all the same, \(I=0.746824132812\ldots\) , and it has to be computed from values of the integrand.

Any such computation uses finitely many samples \(f(x_1),\dots,f(x_n)\) and returns a weighted sum of them. Three questions organize the chapter. Given the sample points, which weights are best, and how large is the error? Given freedom to place the points, where should they go? And, the practical question, how can a program that knows nothing about \(f\) except the values it samples decide that its answer is accurate enough to stop? The last question is the first appearance in this book of a posteriori error control, in which the computation estimates its own error and refines itself where the estimate is large. It works remarkably well, and it can be fooled, and one needs to understand both facts.

6.1 Replacing the integrand by an interpolant

The natural way to obtain a rule is to interpolate \(f\) at the nodes by a polynomial and integrate the polynomial. With \(h=b-a\) and \(m=\tfrac12(a+b)\) , the three simplest rules are

the midpoint rule (integrate the constant interpolant at \(m\) ), the trapezoid rule (the linear interpolant at \(a\) and \(b\) ), and Simpson’s rule (the quadratic interpolant at \(a\) , \(m\) , \(b\) ). The midpoint and trapezoid rules are exact for polynomials of degree \(1\) . Simpson’s rule is exact for quadratics by construction and, by symmetry, also for \((x-m)^3\) , whose integral and whose Simpson sum both vanish; so its degree of exactness is \(3\) , one more than its construction promised.

Errors follow from the interpolation error of Chapter 4. One fact about integrals is needed, a consequence of the intermediate value theorem proved in Real Analysis: if \(g\) is continuous on \([a,b]\) and \(w\) is integrable and does not change sign there, then \(\int_a^bg\,w=g(\xi)\int_a^bw\) for some \(\xi\in[a,b]\) .

Two things stand out. The midpoint rule, with one evaluation, has half the error constant of the trapezoid rule, with two. And the two errors have opposite signs: for a convex integrand the midpoint rule underestimates and the trapezoid rule overestimates, so together they bracket the integral. If the trapezoid error is about \(-2\) times the midpoint error, the combination \(\tfrac13(2M+T)\) should cancel the leading error, and that combination is exactly Simpson’s rule.

6.2 Composite rules and observed order

One panel is rarely enough. Divide \([a,b]\) into \(n\) panels of width \(h=(b-a)/n\) with points \(x_j=a+jh\) and apply a simple rule on each. The composite trapezoid rule is

and for even \(n\) the composite Simpson rule, applied to the \(n/2\) double panels \([x_{2j},x_{2j+2}]\) , is

The exponent drops by one in passing from a single panel to the whole interval, because the number of panels grows like \(1/h\) . The composite trapezoid rule is second-order accurate and the composite Simpson rule fourth-order.

There is an important difference from the difference quotients of Section 1.6. A quadrature sum has positive weights that add up to \(b-a\) , so errors of size \(\epsilon_f\) in the function values change it by at most \((b-a)\epsilon_f\) , however small \(h\) is. Nothing is divided by \(h\) . Numerical integration is a well-conditioned process and numerical differentiation is not, and refining the grid of a quadrature rule does not eventually make things worse.

6.3 Extrapolation: getting more from the trapezoid values

The trapezoid error has more structure than Theorem 6.4 shows. The Euler–Maclaurin formula, which we state without proof, says that for \(f\) with \(2k+2\) continuous derivatives

with universal constants \(\beta_1=\tfrac1{12}\) , \(\beta_2=-\tfrac1{720}\) , and so on. The coefficients do not depend on \(h\) . For the bell-curve integral \(f'(1)-f'(0)=-2e^{-1}\) , and the predicted leading error \(I-T_n\approx\tfrac{1}{12}\cdot2e^{-1}h^2=0.0613\,h^2\) is \(9.58\times10^{-4}\) at \(n=8\) , against the observed \(9.585\times10^{-4}\) .

An expansion in known powers of \(h\) with unknown but fixed coefficients can be exploited. Compute \(T_n\) and \(T_{2n}\) . Their errors are \(c_1h^2+\cdots\) and \(\tfrac14c_1h^2+\cdots\) , so the combination

has no \(h^2\) term at all: its error is of order \(h^4\) . This is Richardson extrapolation, and the combination is in fact the composite Simpson rule \(S_{2n}\) (Exercise 6.7). Repeating the idea on the \(h^4\) terms, then the \(h^6\) terms, gives Romberg integration: with \(R_{k,0}=T_{2^k}\) ,

Each new column gains two orders, and every evaluation of \(f\) is reused, since \(T_{2n}\) needs only the midpoints of the panels of \(T_n\) .

6.4 Gaussian quadrature: choosing the nodes too

An interpolatory rule on \(n\) given nodes is exact for polynomials of degree \(n-1\) . That used the \(n\) weights. If the nodes are also free there are \(2n\) parameters, which suggests that degree \(2n-1\) might be attainable. It is, and the nodes that achieve it are the zeros of an orthogonal polynomial.

The orthogonal polynomials for \([-1,1]\) are the Legendre polynomials: \(P_1=x\) , \(P_2=\tfrac12(3x^2-1)\) , \(P_3=\tfrac12(5x^3-3x)\) . The two-point rule has nodes \(\pm1/\sqrt3\) with weights \(1\) ; the three-point rule has nodes \(0\) and \(\pm\sqrt{3/5}\) with weights \(\tfrac89\) and \(\tfrac59\) . An integral over \([a,b]\) is mapped to \([-1,1]\) by \(x=m+\tfrac12(b-a)t\) , which multiplies the weights by \(\tfrac12(b-a)\) . The error of the two-point rule on \([a,b]\) is \(+\tfrac{h^5}{4320}f''''(\xi)\) , two thirds the size of Simpson’s with two evaluations in place of three.

6.5 Adaptive quadrature: estimating the error as you go

A fixed grid is wasteful when the integrand is smooth in most places and difficult in a few. An adaptive method decides locally. On a panel \([a,b]\) compute Simpson’s rule once, \(S_1\) , and then on its two halves, \(S_2\) . If \(f''''\) varies little across the panel, Theorem 6.3 gives \(I-S_1\approx-\tfrac{h^5}{2880}f''''\) and \(I-S_2\approx-2\cdot\tfrac{(h/2)^5}{2880}f''''=\tfrac1{16}(I-S_1)\) . Subtracting, \(S_2-S_1\approx\tfrac{15}{16}(I-S_1)=15\,(I-S_2)\) , so

The difference of two computed numbers estimates the error of the better one. The algorithm accepts the panel if \(\lvert S_2-S_1\rvert\le15\,\tau\) , where \(\tau\) is the tolerance allotted to that panel, and otherwise bisects the panel and applies itself to each half with tolerance \(\tau/2\) . On acceptance it may as well add the estimated error to \(S_2\) .

def adaptive_simpson(f, a, b, tol):
    def refine(a, fa, m, fm, b, fb, whole, tol):
        lm, rm = (a + m) / 2, (m + b) / 2
        flm, frm = f(lm), f(rm)
        left = (m - a) / 6 * (fa + 4 * flm + fm)
        right = (b - m) / 6 * (fm + 4 * frm + fb)
        diff = left + right - whole
        if abs(diff) <= 15 * tol:
            return left + right + diff / 15
        return (refine(a, fa, lm, flm, m, fm, left, tol / 2)
                + refine(m, fm, rm, frm, b, fb, right, tol / 2))
    fa, fb, m = f(a), f(b), (a + b) / 2
    fm = f(m)
    return refine(a, fa, m, fm, b, fb, (b - a) / 6 * (fa + 4 * fm + fb), tol)

Each call costs two new evaluations of \(f\) , and every value is passed down so that nothing is computed twice. A production version adds a limit on the depth of recursion.

Example 6.6: A square-root singularity

The integral \(\int_0^1\sqrt x\,dx=\tfrac23\) has an integrand whose derivatives are unbounded at \(0\) . Theorem 6.4 does not apply, and the composite Simpson rule converges only like \(h^{3/2}\) : its errors are \(2.48\times10^{-6}\) with \(1024\) panels, \(3.10\times10^{-7}\) with \(4096\) , and \(3.87\times10^{-8}\) with \(16384\) , a factor of \(8\) for each fourfold increase in work.

The adaptive routine with tolerance \(10^{-6}\) returns \(0.666666647130\) , with error \(1.95\times10^{-8}\) , using \(105\) function evaluations. Of those \(105\) points, \(62\) lie in \([0,0.01]\) , and the closest pair of samples is \(2^{-23}\approx1.2\times10^{-7}\) apart; on the right half of the interval, where \(\sqrt x\) is smooth, the accepted panels are a quarter of the interval wide. To reach the same accuracy on a uniform grid would take more than \(16000\) evaluations. The adaptive method found the singularity without being told where it was.

Figure 6.1 summarizes the two routes to efficiency in this chapter. On the left the integrand is \(1/((x-0.3)^2+0.001)+1/((x-0.9)^2+0.04)-6\) , with a sharp peak at \(0.3\) and a broad one at \(0.9\) , and the adaptive samples crowd around the sharp peak while a handful suffice for the rest. On the right the integrand is smooth, the slopes of the two composite rules display their orders exactly as in Figure 1.1, and the Gauss rules converge faster than any fixed power.

6.6 When the error estimate lies

The estimate \((S_2-S_1)/15\) rests on an assumption, that the fourth derivative is nearly constant across the panel. The program cannot verify the assumption. It sees five numbers.

Exercises

1. Compute the midpoint, trapezoid, and Simpson approximations to \(\int_1^2dx/x=\ln2\) by hand, find the three errors, and check their signs and sizes against Theorems 6.2 and 6.3.

   Solution

   With \(f(x)=1/x\) , \(h=1\) , and \(m=1.5\) : \(M=f(1.5)=0.666667\) , \(T=\tfrac12(1+\tfrac12)=0.75\) , and \(S=\tfrac16(1+4/1.5+0.5)=0.694444\) . Since \(\ln2=0.693147\) , the errors \(I-Q\) are \(0.02648\) for the midpoint rule, \(-0.05685\) for the trapezoid rule, and \(-0.001297\) for Simpson’s rule.

   Signs. The integrand is convex, with \(f''(x)=2/x^3>0\) . Theorem 6.2 gives \(I-M=\tfrac1{24}f''(\eta)>0\) and \(I-T=-\tfrac1{12}f''(\xi)<0\) : the midpoint rule underestimates and the trapezoid rule overestimates, as observed. Also \(f''''(x)=24/x^5>0\) , so Theorem 6.3 gives \(I-S<0\) .

   Sizes. On \([1,2]\) , \(f''\) ranges over \([0.25,2]\) , so \(I-M\) must lie in \([0.0104,0.0833]\) and \(I-T\) in \([-0.1667,-0.0208]\) ; both observed errors do. And \(f''''\) ranges over \([0.75,24]\) , so \(I-S\) must lie in \([-0.00833,-0.00026]\) , which contains \(-0.001297\) . The trapezoid error is \(-2.15\) times the midpoint error, near the ratio \(-2\) of the leading terms; the two would be exactly in that ratio if \(f''\) were constant.

2. Verify directly that Simpson’s rule on \([-1,1]\) integrates \(x^3\) exactly but not \(x^4\) , compute its error for \(x^4\) , and confirm that the error agrees with Theorem 6.3.

   Answer

   On \([-1,1]\) , \(S(f)=\tfrac13\bigl(f(-1)+4f(0)+f(1)\bigr)\) . For \(x^3\) : \(S=\tfrac13(-1+0+1)=0=\int_{-1}^{1}x^3dx\) . For \(x^4\) : \(S=\tfrac23\) while \(\int_{-1}^{1}x^4dx=\tfrac25\) , so the error is \(\tfrac25-\tfrac23=-\tfrac4{15}\) . Theorem 6.3 with \(h=2\) and \(f''''=24\) gives \(-\tfrac{32}{2880}\cdot24=-\tfrac4{15}\) , in agreement.

3. Find weights \(w_0,w_1,w_2\) such that \(\int_0^1f(x)\,dx\approx w_0f(0)+w_1f(\tfrac13)+w_2f(1)\) is exact for all quadratics, and determine the degree of exactness of the resulting rule.

   Solution

   Require exactness for \(1\) , \(x\) , and \(x^2\) , whose integrals over \([0,1]\) are \(1\) , \(\tfrac12\) , and \(\tfrac13\) :

   \[ w_0+w_1+w_2=1,\qquad\tfrac13w_1+w_2=\tfrac12,\qquad\tfrac19w_1+w_2=\tfrac13 . \]

   Subtracting the third equation from the second gives \(\tfrac29w_1=\tfrac16\) , so \(w_1=\tfrac34\) ; then \(w_2=\tfrac12-\tfrac14=\tfrac14\) and \(w_0=0\) . The rule is

   \[ \int_0^1f(x)\,dx\approx\tfrac34f(\tfrac13)+\tfrac14f(1). \]

   By linearity it is exact for every quadratic. For \(x^3\) it gives \(\tfrac34\cdot\tfrac1{27}+\tfrac14=\tfrac5{18}\) , whereas the integral is \(\tfrac14\) , so the degree of exactness is \(2\) .

   The node \(0\) dropped out: its weight is zero. What remains is a two-point rule with degree of exactness \(2\) , one more than a two-point rule on arbitrary nodes can be expected to have, because the node \(\tfrac13\) happens to be the right partner for a node fixed at the endpoint \(1\) . Rules of this kind, with one node pinned at an endpoint and the others placed optimally, are called Gauss–Radau rules; they sit between the interpolatory rules on given nodes and the Gauss rules of Theorem 6.5.

4. A rule \(Q(f)=\sum_{i=1}^{n}w_if(x_i)\) with distinct nodes is exact for all polynomials of degree at most \(n-1\) . Prove that its weights must be \(w_i=\int_a^b\ell_i(x)\,dx\) , where \(\ell_i\) are the Lagrange basis polynomials for the nodes, so that the rule is the integral of the interpolating polynomial.

   Solution

   Each Lagrange basis polynomial \(\ell_j\) has degree \(n-1\) , so by hypothesis the rule integrates it exactly:

   \[ \int_a^b\ell_j(x)\,dx=Q(\ell_j)=\sum_{i=1}^{n}w_i\,\ell_j(x_i)=w_j , \]

   because \(\ell_j(x_i)\) equals \(1\) when \(i=j\) and \(0\) otherwise. That determines every weight.

   Now let \(f\) be any function defined at the nodes and let \(p=\sum_if(x_i)\ell_i\) be its interpolating polynomial of degree at most \(n-1\) . Then

   \[ Q(f)=\sum_iw_if(x_i)=\sum_if(x_i)\int_a^b\ell_i(x)\,dx=\int_a^bp(x)\,dx . \]

   So a rule on \(n\) distinct nodes is exact for polynomials of degree at most \(n-1\) if and only if it is obtained by integrating the interpolant: the converse direction holds because the interpolant of a polynomial of degree at most \(n-1\) is that polynomial itself. In particular the nodes determine the weights, and the only way to gain further degrees of exactness is through the placement of the nodes, which is the idea of Gaussian quadrature.

5. Use Theorem 6.4 to find the number of panels that guarantees an error below \(10^{-6}\) for \(\int_0^\pi\sin x\,dx\) with the composite trapezoid rule and with the composite Simpson rule.

   Answer

   Trapezoid: \(\tfrac{\pi h^2}{12}\cdot1\le10^{-6}\) gives \(h\le1.954\times10^{-3}\) , so \(n\ge1607.4\) : use \(1608\) panels. Simpson: \(\tfrac{\pi h^4}{180}\le10^{-6}\) gives \(h\le0.0870\) , so \(n\ge36.1\) ; since \(n\) must be even, use \(38\) .

6. Without using the exact value of the integral, estimate the order of the composite trapezoid rule from the three values \(T_2\) , \(T_4\) , \(T_8\) quoted in Example 6.2, and explain why the estimate is legitimate.

   Answer

   \((T_2-T_4)/(T_4-T_8)=0.0116138/0.0028815=4.03\) , so \(p=\log_24.03=2.01\) . If \(T(h)=I+Ch^{p}+\cdots\) , then \(T(h)-T(h/2)=Ch^{p}(1-2^{-p})+\cdots\) , and the unknown \(I\) cancels from each difference while the unknown \(C\) cancels from their ratio, leaving \(2^{p}\) . The estimate is legitimate once \(h\) is small enough for the leading term to dominate, which is tested by computing the ratio again at the next level and seeing it stay put.

7. Prove that \(\tfrac13\bigl(4T_{2n}(f)-T_n(f)\bigr)=S_{2n}(f)\) for every function \(f\) and every \(n\) .

   Solution

   Let \(h\) be the panel width of \(T_{2n}\) , so that \(T_n\) has panel width \(2h\) and uses the even-numbered points \(x_0,x_2,\dots,x_{2n}\) . Writing \(f_j=f(x_j)\) ,

   \[ T_{2n}=h\bigl(\tfrac12f_0+f_1+f_2+\dots+f_{2n-1}+\tfrac12f_{2n}\bigr),\qquad T_n=2h\bigl(\tfrac12f_0+f_2+f_4+\dots+f_{2n-2}+\tfrac12f_{2n}\bigr). \]

   Form \(4T_{2n}-T_n\) and collect the coefficient of each \(f_j\) . The endpoints get \(4\cdot\tfrac12h-2h\cdot\tfrac12=h\) . An odd-numbered point appears only in \(T_{2n}\) and gets \(4h\) . An even-numbered interior point gets \(4h-2h=2h\) . Therefore

   \[ \frac{4T_{2n}-T_n}{3}=\frac{h}{3}\bigl(f_0+4f_1+2f_2+4f_3+\dots+2f_{2n-2}+4f_{2n-1}+f_{2n}\bigr)=S_{2n}(f). \]

   The identity is purely algebraic and holds for every \(f\) ; no smoothness is used. It shows that Simpson’s rule is the first Richardson extrapolation of the trapezoid rule, and it explains why Simpson’s rule gains two orders: the extrapolation removes the \(h^2\) term of the Euler–Maclaurin expansion.

8. Build the Romberg table from \(T_1\) , \(T_2\) , \(T_4\) for \(\int_0^14/(1+x^2)\,dx=\pi\) , giving each entry to six decimals and the errors of the entries in the last row.

   Answer

   \(T_1=3.000000\) , \(T_2=3.100000\) , \(T_4=3.131176\) . First extrapolation: \(R_{1,1}=3.133333\) and \(R_{2,1}=3.141569\) . Second: \(R_{2,2}=(16R_{2,1}-R_{1,1})/15=3.142118\) . The errors \(\pi-R\) in the last row are \(1.04\times10^{-2}\) , \(2.4\times10^{-5}\) , and \(-5.2\times10^{-4}\) . The second extrapolation made things worse. For this integrand \(f'''(1)=f'''(0)=0\) , so the \(h^4\) term of the Euler–Maclaurin expansion vanishes, Simpson’s rule is already better than fourth order, and removing a term that is not there only injects error from the coarser entry.

9. Determine the two-point Gauss rule on \([-1,1]\) by requiring a symmetric rule \(f(-\alpha)+f(\alpha)\) to be exact for \(x^2\) , verify its degree of exactness, and compare its error with that of Simpson’s rule for \(\int_{-1}^{1}\cos x\,dx\) .

   Solution

   A symmetric rule \(Q(f)=f(-\alpha)+f(\alpha)\) integrates every odd power exactly, since both sides vanish, and integrates constants exactly, since \(1+1=2\) . Exactness for \(x^2\) requires \(2\alpha^2=\tfrac23\) , so \(\alpha=1/\sqrt3\) . The rule is then exact for \(1,x,x^2,x^3\) , hence for all cubics. For \(x^4\) it gives \(2\alpha^4=\tfrac29\) , not \(\tfrac25\) , so the degree of exactness is \(3\) , which is \(2n-1\) for \(n=2\) , as Theorem 6.5 asserts.

   For \(\int_{-1}^{1}\cos x\,dx=2\sin1=1.682942\) : the Gauss rule gives \(2\cos(1/\sqrt3)=1.675824\) , with error \(7.12\times10^{-3}\) . Simpson’s rule gives \(\tfrac13(\cos1+4+\cos1)=1.693535\) , with error \(-1.059\times10^{-2}\) . The Gauss rule, with two evaluations, is more accurate than Simpson’s rule with three, and the errors have opposite signs. Their ratio in magnitude is \(1.49\) , matching the ratio \(4320/2880=1.5\) of the constants in the two error formulas, both of which involve \(f''''=\cos\) at an intermediate point.

10. Map the three-point Gauss–Legendre rule to \([0,1]\) , giving nodes and weights, and use it to approximate \(\int_0^1dx/(1+x)=\ln2\) ; report the error.

    Answer

    The map \(x=\tfrac12+\tfrac12t\) sends the nodes \(0,\pm\sqrt{3/5}\) to \(0.112702\) , \(0.5\) , \(0.887298\) and halves the weights to \(\tfrac5{18},\tfrac49,\tfrac5{18}\) . The rule gives \(0.693122\) for \(\int_0^1dx/(1+x)\) , an error of \(2.5\times10^{-5}\) against \(\ln2=0.693147\) .

11. Prove that if \(f\) is convex and twice continuously differentiable on \([a,b]\) then \(M(f)\le\int_a^bf\le T(f)\) , and state the corresponding result for the composite rules.

    Hint

    Convexity with two continuous derivatives means \(f''\ge0\) . Read the signs of the two errors in Theorem 6.2. For the composite rules, add the single-panel inequalities over the panels to get \(M_n(f)\le\int_a^bf\le T_n(f)\) , a computable two-sided bound.

12. Doubling the number of panels of the composite trapezoid rule for \(\int_0^1\sqrt x\,dx\) reduces the error by a factor of about \(2.8\) instead of \(4\) . Identify the hypothesis of Theorem 6.4 that fails, say what order the factor indicates, and give a change of variable after which the rule converges at its normal rate.

    Answer

    Theorem 6.4 requires \(f''\) to be continuous on \([0,1]\) ; for \(\sqrt x\) it is unbounded at \(0\) . A factor of \(2.8=2^{1.5}\) per doubling indicates order \(\tfrac32\) . With \(x=t^2\) the integral becomes \(\int_0^12t^2\,dt\) , whose integrand is a polynomial, and the trapezoid rule converges at order \(2\) .

13. Use the Euler–Maclaurin formula to explain the behaviour in Example 6.4, and predict what happens to the trapezoid rule for \(\int_0^{2\pi}e^{\cos x}dx\) if the interval is changed to \([0,\pi]\) .

    Hint

    For a smooth integrand of period \(2\pi\) integrated over a full period, every odd derivative takes the same value at both ends, so every coefficient \(c_j\) in the expansion vanishes. On \([0,\pi]\) the integrand \(e^{\cos x}\) is still even about both endpoints, so its odd derivatives vanish there; check whether that again kills every \(c_j\) .

14. Derive the estimate \(I-S_2\approx(S_2-S_1)/15\) of Section 6.5 from Theorem 6.3, stating the assumption on \(f''''\) explicitly, and derive the analogous estimate for the trapezoid rule on one panel and its two halves.

    Solution

    Simpson. Let the panel have width \(h\) and assume \(f''''\) is approximately a constant \(F\) on it. Theorem 6.3 gives \(I-S_1\approx-\tfrac{h^5}{2880}F\) . The refined value \(S_2\) is the sum of two Simpson rules on half-panels of width \(h/2\) , each with error approximately \(-\tfrac{(h/2)^5}{2880}F\) , so

    \[ I-S_2\approx-2\cdot\frac{h^5}{32\cdot2880}F=\tfrac1{16}(I-S_1). \]

    Subtracting, \(S_2-S_1=(I-S_1)-(I-S_2)\approx16(I-S_2)-(I-S_2)=15(I-S_2)\) , hence \(I-S_2\approx(S_2-S_1)/15\) .

    Trapezoid. Assume \(f''\approx G\) on the panel. Theorem 6.2 gives \(I-T_1\approx-\tfrac{h^3}{12}G\) and, for two half-panels, \(I-T_2\approx-2\cdot\tfrac{(h/2)^3}{12}G=\tfrac14(I-T_1)\) . Then \(T_2-T_1\approx4(I-T_2)-(I-T_2)=3(I-T_2)\) , so \(I-T_2\approx(T_2-T_1)/3\) .

    In both cases the assumption is that the relevant derivative does not change much across the panel, so that the unknown points \(\xi\) in the error formulas for the whole panel and for its halves may be treated as giving the same value. When that fails, as in Example 6.7, the estimate fails with it. The divisors \(15=2^4-1\) and \(3=2^2-1\) are the Richardson factors for fourth and second order.

15. For \(f(x)=x^4\) on \([0,1]\) compute \(S_1\) on the whole interval and \(S_2\) on its two halves by hand, form the estimate \((S_2-S_1)/15\) , and compare it with the true error of \(S_2\) ; explain why the agreement is exact.

    Solution

    On \([0,1]\) : \(S_1=\tfrac16\bigl(0+4\cdot\tfrac1{16}+1\bigr)=\tfrac5{24}=0.208333\) . On the halves, each of width \(\tfrac12\) :

    \[ \tfrac1{12}\bigl(0+4\cdot\tfrac1{256}+\tfrac1{16}\bigr)=\tfrac{5}{768},\qquad\tfrac1{12}\bigl(\tfrac1{16}+4\cdot\tfrac{81}{256}+1\bigr)=\tfrac{149}{768}, \]

    so \(S_2=\tfrac{154}{768}=\tfrac{77}{384}=0.200521\) . The exact integral is \(\tfrac15\) . The estimate is

    \[ \frac{S_2-S_1}{15}=\frac1{15}\Bigl(\frac{77}{384}-\frac{80}{384}\Bigr)=-\tfrac1{1920}=-5.208\times10^{-4}, \]

    and the true error is \(I-S_2=\tfrac15-\tfrac{77}{384}=-\tfrac1{1920}\) . They agree exactly.

    The reason is that for \(f(x)=x^4\) the fourth derivative is the constant \(24\) , so the assumption behind the estimate holds with no approximation: Theorem 6.3 gives \(I-S_1=-\tfrac{24}{2880}=-\tfrac1{120}\) exactly and \(I-S_2=-2\cdot\tfrac{24}{32\cdot2880}=-\tfrac1{1920}\) exactly, and the derivation of Exercise 6.14 becomes an identity. Adding the estimate to \(S_2\) therefore reproduces the integral exactly; in general the corrected value \(S_2+(S_2-S_1)/15\) is exact for all polynomials of degree at most \(5\) .

16. An integrand equals \(1\) except on an interval of width \(10^{-3}\) centred at \(0.37\) , where it rises smoothly to \(1000\) . Explain why adaptive Simpson started on \([0,1]\) with a loose tolerance is likely to return \(1\) , why no method based on sampling can guarantee to detect the spike, and what a user who knows about the spike should do.

    Hint

    Locate the first five sample points and ask whether any lies within \(5\times10^{-4}\) of \(0.37\) . If none does, what are \(S_1\) and \(S_2\) ? For the second part, given any finite set of sample points one can place a spike between them. A user who knows the spike is there should split the integral at \(0.3695\) and \(0.3705\) and integrate the three pieces separately.

17. Let the computed values \(\hat f(x_i)\) satisfy \(\lvert\hat f(x_i)-f(x_i)\rvert\le\epsilon_f\) . Show that the composite trapezoid and Simpson sums change by at most \((b-a)\epsilon_f\) , and contrast this with the rounding term of Theorem 1.11 to explain why quadrature has no optimal step size.

    Solution

    Let \(Q(f)=\sum_iw_if(x_i)\) be either composite rule. Its weights are positive: \(\tfrac h2\) and \(h\) for the trapezoid rule, and \(\tfrac h3\) , \(\tfrac{4h}3\) , \(\tfrac{2h}3\) for Simpson’s rule. Both rules integrate the constant function \(1\) exactly, so \(\sum_iw_i=b-a\) . Therefore

    \[ \bigl\lvert Q(\hat f)-Q(f)\bigr\rvert=\Bigl\lvert\sum_iw_i\bigl(\hat f(x_i)-f(x_i)\bigr)\Bigr\rvert\le\epsilon_f\sum_i\lvert w_i\rvert=\epsilon_f\sum_iw_i=(b-a)\,\epsilon_f . \]

    The bound does not depend on \(h\) or on the number of panels. Combined with Theorem 6.4, the total error of the computed composite rule is at most \(Ch^{p}+(b-a)\epsilon_f\) , a quantity that decreases as \(h\) decreases and levels off at \((b-a)\epsilon_f\) . There is no step below which refinement makes the answer worse, only a floor beyond which it stops helping.

    For the difference quotient of Theorem 1.11 the corresponding term is \(2\epsilon_f/h\) , which grows without bound as \(h\to0\) : the weights of a differentiation formula are \(\pm1/h\) , of mixed sign and unbounded sum of magnitudes. The contrast is one of conditioning. Integration averages errors in the data; differentiation amplifies them. Positive weights are the reason, and they are also why Gauss rules, whose weights are positive by Theorem 6.5, are as benign as the composite rules.

18. The integral \(\int_0^1\cos x/\sqrt x\,dx\) has an integrable singularity at \(0\) . Use the substitution \(x=t^2\) to rewrite it with a smooth integrand, apply the composite Simpson rule with \(n=4\) and \(n=8\) , and estimate the error of the second value from the two.

    Solution

    With \(x=t^2\) , \(dx=2t\,dt\) and \(\sqrt x=t\) , so

    \[ \int_0^1\frac{\cos x}{\sqrt x}\,dx=\int_0^12\cos(t^2)\,dt , \]

    and the new integrand \(g(t)=2\cos(t^2)\) has derivatives of all orders on \([0,1]\) . The composite Simpson rule gives \(S_4=1.8090025\) and \(S_8=1.8090483\) .

    By the Richardson reasoning of Section 6.5 the error of \(S_8\) is estimated by \((S_8-S_4)/15=3.1\times10^{-6}\) . The integral is in fact \(1.8090485\) , so the true error of \(S_8\) is \(1.6\times10^{-7}\) , twenty times smaller than the estimate. The estimate assumed that halving \(h\) divides the error by \(16\) ; here it divided it by nearly \(300\) , because for this integrand the leading error coefficient is small and the error changes sign between \(n=8\) and \(n=16\) , so the computation is not yet in the asymptotic regime. An error estimate can err in either direction. This time it was conservative, which is the harmless direction, but the episode shows why agreement of successive estimates, not a single one, is the evidence to look for.

19. Simpson’s rule with \(n=2\) and the two-point Gauss rule were both applied to \(\int_0^1e^{-x^2}dx\) in this chapter. Show that the ratio and the signs of their errors agree with the two error formulas \(-\tfrac{h^5}{2880}f''''\) and \(+\tfrac{h^5}{4320}f''''\) , and find the weighted combination of the two rules that cancels the leading error.

    Answer

    \(I-S=-3.563\times10^{-4}\) and \(I-G_2=+2.294\times10^{-4}\) : opposite signs, with ratio \(-1.55\) , close to \(-4320/2880=-1.5\) . The combination \(\tfrac15(2S+3G_2)\) has leading error \(\tfrac15\bigl(-\tfrac{2}{2880}+\tfrac{3}{4320}\bigr)h^5f''''=0\) ; numerically it gives \(0.7468290\) , with error \(-4.9\times10^{-6}\) , about seventy times smaller than Simpson’s.

20. Integrating the piecewise linear interpolant of Section 5.1 gives the composite trapezoid rule. Use Proposition 5.1 to derive the bound \((b-a)h^2\lVert f''\rVert_\infty/8\) for its error, compare with the constant \(\tfrac1{12}\) of Theorem 6.4, and explain why integrating a pointwise bound loses a factor.

    Hint

    Integrate the inequality \(\lvert f(x)-s(x)\rvert\le\tfrac18h^2\lVert f''\rVert_\infty\) over \([a,b]\) . The bound assumes the worst-case error at every point simultaneously, whereas the true error on a panel is \(\tfrac12f''(\xi_x)(x-x_j)(x-x_{j+1})\) , which vanishes at both ends of the panel; its average over the panel is two thirds of its maximum, and \(\tfrac23\cdot\tfrac18=\tfrac1{12}\) .