Math 128a - Spring 2002 - Notes on Integration Read sections 7.2 to 7.5 of the text. Integration, or Quadrature We will have several methods for approximating integral (from a to b) f(x)dx, depending on what we know (or don't know) about f(x). In particular, the choice of method depends on how smooth f(x) is (how many derivatives does it have) where we can evaluate f(x) a fixed number of points that are given to us, versus points that we can choose if the points are given to us, are they regularly or irreguarly spaced? We note that this field has lots of names attached to various formulas for historical reasons; I prefer that you remember the ideas behind the formulas instead of the names, if you can't remember everything. Integration via Lagrange Interpolation: Recall the Lagrange interpolation formula for f(x) on x{0},...,x{n} in [a,b]: p(x) = sum (from i=0 to n) f(x{i})*L{i}(x) where L{i}(x) is a Lagrange polynomial (1 at x{i}, 0 at other x{j}). Then integral (from a to b) f(x)dx ~ integral (from a to b) p(x)dx = integral (from a to b) ( sum (from i=0 to n) f(x{i})*L{i}(x)dx ) = sum (from i=0 to n) integral (from a to b) f(x{i})*L{i}(x)dx = sum (from i=0 to n) f(x{i}) * integral (from a to b) L{i}(x)dx = sum (from i=0 to n) f(x{i}) * A{i} where A(i) is a constant, that depends on x{0},...,x{n},a and b, but not f(x). This means that A(i) can be precomputed. If the x{i} are evenly spaced, this is called a Newton-Cotes formula. An error term can be gotten by integrating the error term for interpolation. Ex: Trapezoidal Rule: we take x{0}=a and x{1}=b, fit a straight line (this approximates the area under the curve with the area of a trapezoid, hence the name). One can easily compute the formula: p(x) = (b-a)/2 * (f(a)+f(b)) The error is integral (from a to b) f''(z(t))/2*(x-a)*(x-b) dx = f''(z) * integral (from a to b) (x-a)*(x-b)/2 dx for some z in [a,b] by the Mean Value Theorem for Integrals. = -f''(z)*(b-a)^3/12 Change of interval of integration: It appears that we need new coefficients for each interval of integration [a,b]. We can save ourselves work by changing the variable of integration from from x to t where x = (b-a)*t + a, where t ranges from 0 to 1: integral (from a to b) f(x)dx = integral (from 0 to 1) f((b-a)*t+a) (b-a)*dt ~ (b-a) * sum (from i=0 to n) f((b-a)*t{i}+a)*A{i} where t{i} are in [0,1]. Thus any quadrature formula for [0,1] can easily be turned into a formula for [a,b]. In what follows, we simplify by looking at [0,1]. Method of undetermined coefficients to determine the A{i} If we use n+1 interpolation points, then we know that p(x) will exactly equal any polynomial of degree n or less, and so the integration formula integral (from 0 to 1) f(x)dx ~ sum (from i=0 to n) f(x{i})*A{i} should be exact for f(x) any polynomial of degree n or less. This means integral (from 0 to 1) x^j dx = 1/(j+1) = sum (from i=0 to n) x{i}^j * A{i} for j=0 to n. This is an n+1-by-n+1 Vandermonde systems of linear equations Ex: Simpson's rule: This corresponds to using x{0}=0, x{1}=1/2 and x{2}=1. Applying above method yields integral (0 to 1) f(x)dx = (1/6)*(f(0) + 4*f(1/2) + f(1)), Here is one way to figure out the error. Assume we are integrating from a to b=a+2*h, so the formula is integral (a to a+2h) f(x) = (h/3)*(f(a) + 4*f(a+h) + f(a+2*h)) First, do a Taylor expansion of the right hand side to get (h/3)*(f(a) + 4*f(a+h) + f(a+2*h) = (h/3)*(f(a) + 4*(f(a) + h*f'(a) + h^2/2*f''(a) + h^3/6*f'''(a) + h^4/24*f^(4)(a) + O(h^5) ) + f(a) + 2*h*f'(a) + 4*h^2/2*f''(a) + 8*h^3/6*f'''(a) + 16*h^4/24*f^(4)(a) + O(h^5) ) = 2*h*f(a) + 2*h^2*f'(a) + (4/3)*h^3*f''(a) + (2/3)*h^4*f'''(a) + (5/18)*h^5*f^(4)(a) + O(h^6) Second, in the left hand side integrate the Taylor expansion of f(x) to get integral (a to a+2h) f(x)dx = 2*h*f(a) + 2*h^2*f'(a) + (4/3)*h^3*f''(a) + + (2/3)*h^4*f'''(a) + (4/15)*h^5*f^(4)(a) + O(h^6) Comparing, we get integral (a to a+2h) f(x)dx = (h/3)*(f(a) + 4*f(a+h) + f(a+2*h) + (4/15 - 5/18)*h^5*f^(4)(a) + O(h^6) = (h/3)*(f(a) + 4*f(a+h) + f(a+2*h) - (1/90)*h^5*f^(4)(a) + O(h^6) Summarizing, the error in Simpson's rule when integrating from a to b is -h^5/90*f^(4)(a) = -((b-a)/2)^5/90*f^(4)(z) = -(b-a)^5/2880*f^(4)(z) Composite Integration (via piecewise Lagrange Interpolation): This refers to dividing the interval [a,b] into equal subintervals, integrating via Lagrange interpolation on each subinterval independently, and adding the results. Ex: Composite Trapezoial Rule: We sample [0,1] at x{i}=i*h, h=1/n, so x{n}=1, and use the Trapezoidal rule on each interval [x{i},x{i+1}]. This yields integral (from 0 to 1) f(x)dx = sum (j=0 to n-1) integral (from j*h to (j+1)*h) f(x)dx ~ sum (j=0 to n-1) (h/2)*(f(j*h)+f(j+1)*h) = h * sum'' (j=0 to n) f(j*h) where sum'' means that the j=0 and j=n terms are multiplied by 1/2. The error is the sum of the errors, and can be simplified: sum (j=0 to n-1) -f''(z{j})*h^3/12 = -h^2/12 * sum (j=0 to n-1) f''(z{j})*h = -h^2/12 * f''(z) where the last simplification is due to the Intermediate Value Theorem Ex: Composite Simpson: The same approach, if we have an even number of intervals (since the basic Simpson's rule requires 2 intervals) is: integral (from 0 to 1) f(x)dx = sum (j=0 to n-2 step 2) integral (from j*h to (j+2)*h) f(x)dx ~ sum (j=0 to n-2 step 2) h/3*f(j*h)+4*h/3*f(j+1)*h)+h/3*f((j+2)*h) = (h/3) * [ f(0) + f(1) ] + (2*h/3) * sum (j=2 to n-2 step 2) f(j*h) + (4*h/3) * sum (j=0 to n-2 step 2) f((j+1)*h) + The same approach to summing the errors for each interval yields -h^4/180*f^(4)(z) Adaptive Quadrature: To build a reliable "black box" integrator, one should (1) let the user specify any desired accuracy tol: | true integral of f - computed integral of f | <= tol (2) accommodate functions f(x) that are not very smooth. To achieve (1), one needs to compute an error estimate as one proceeds, and recompute the integral if the error estimate is too large. To achieve (2), one cannot use a method that assumes f(x) has derivatives of high order. But we have seen that one still needs some derivatives to get any error estimate at all. Our compromise will be to use Simpson's rule, whose error bound assumes a fourth derivative: integral (a to b) f(x)dx = ((b-a)/6)*(f(a) + 4*f((a+b)/2) + f(b)), -(b-a)^5/2880*f^(4)(z) = S(a,b) + E(a,b) Here S(a,b) = ((b-a)/6)*(f(a) + 4*f((a+b)/2) + f(b)) is the formula for Simpson's rule and E(a,b) = -(b-a)^5/2880*f^(4)(z) is the error. Actually, we will use a composite Simpson's rule on an adaptively chosen set of points a=x{0} < x{1} < ... < x{n} = b to get integral (from a to b) f(x) dx = sum (from i=0 to n-1) S(x{i},x{i+1}) + e(x{i},x{i+1}) and choose the x{i} adaptively so that | sum (from i=0 to n-1) e(x{i},x{i+1}) | <= tol Here is the overall algorithm. We state it first recursively, and then using an explicit data structure called a Stack to keep track of which intervals we are working on. Just think of the Stack like a stack of dishes, where you can put new dishes on top, and when you want to remove something, you get the last dish you put on top of the stack. Assume that we have an estimator e(a,b) ~ E(a,b) (we will describe it below) Adaptive Quaduature (Recursive) V = 0, TotalError = 0 [V,TotalError] = AdaptQuad(a,b,S(a,b),e(a,b),V,TotalError) ... note: could make V and TotalError global variables function [V,TotalError] = AdaptQuad(low,up,S,e,V,TotalError) if e not small enough mid = (low+up)/2 [V,TotalError] = AdaptQuad(low,mid,S(low,mid),e(low,mid),V,TotalError) [V,TotalError] = AdaptQuad(mid,up,S(mid,up),e(mid,up),V,TotalError) else V = V + S TotalError = TotalError + e endif Adaptive Quaduature (Nonrecursive) V = 0; TotalError = 0; Put [a,b,S(a,b),e(a,b)] on the Stack. while there are items on the Stack remove [low,up,S,e] from the Stack if e not small enough ... split the interval [low,up] in half mid = (low+up)/2; put [low,mid,S(low,mid),e(low,mid)] on the Stack put [mid,up,S(mid,up),e(mid,up)] on the Stack else V = V + S TotalError = TotalError + e end end while To make this work, we need to describe how to estimate the error e(a,b), and how to decide whether it is small enough. First, we describe how to decide if it is small enough. Suppose that the final list of points determined by the algorithm is a=x{0},x{1},...,x{n}=b. Then if e(x{i},x{i+1}) <= tol*(x{i}-x{i+1})/(b-a) then the TotalError will be at most sum (from i=0 to n-1) e(x{i},x{i+1}) <= sum (from i=0 to n-1) tol*(x{i}-x{i+1})/(b-a) = tol*(b-a)/(b-a) = tol as desired. So our test "if e not small enough" becomes "if ( |e| >= tol*(up-low)/(b-a) )" Here is how we estimate e(a,b): we estimate integral (a to b) f(x)dx two ways (here c = (a+b)/2): (*) integral (a to b) f(x)dx = S(a,b) + E(a,b) = S(a,b) -(b-a)^5/2880*f^(4)(z1) and (**) integral (a to b) f(x)dx = integral (a to c) f(x)dx + integral (c to b) f(x)dx = S(a,c) + E(a,c) + S(c,b) + E(c,b) = S(a,c) + S(c,b) -(c-a)^5/2880*f^(4)(z2) -(b-c)^5/2880*f^(4)(z3) = S(a,c) + S(c,b) -(b-a)^5/2^5/2880*f^(4)(z2) -(b-a)^5/2^5/2880*f^(4)(z3) = S(a,c) + S(c,b) -(b-a)^5/16/2880 * [ f^(4)(z2) + f^(4)(z3) ] /2 = S(a,c) + S(c,b) -(b-a)^5/16/2880 * f^(4)(z4) by the Intermediate Value Theorem Now we don't know much about f^(4), but we make the assumption that f(x) is smooth enough it changes slowly, so that f^(4)(z4)~f^(4)(z1), since z1 and z4 are close together. Now we can set (*) = (**) and get S(a,b) -(b-a)^5/2880*f^(4)(z1) = S(a,c) + S(c,b) -(b-a)^5/16/2880 * f^(4)(z4) or S(a,b) + E = S(a,c) + S(c,b) + E/16 or E = ( S(a,c) + S(c,b) - S(a,b) ) * 16/15 This is our error estimator e(a,b). This trick is essentially the same trick as we used in Richard Extrapolation: we solved the problem with an interval of width h=b-a, and 2 intervals of length h/2, and combined the estimates to solve for the error. This is now enough to fill in the details of the Adaptive Quaduature algorithm above. There is one obvious improvement to make: Computing an error bound e(a,b) required us to compute the solution a little more accurately, namely as S(a,c)+S(c,b). In fact we have an error bound for it too, e(a,b)/16, so we expect it to be 16 times more accurate. Therefore, we should use S(a,c)+S(c,b) as the final estimate of integral (a to b) f(x)dx and ( S(a,c) + S(c,b) - S(a,b) ) /15 as the error bound e(a,b). So here is the final Adaptive Quadrature algorithm. It uses the function S(a,b) = ((b-a)/6)*(f(a) + 4*f((a+b)/2) + f(b)) V = 0 TotalError = 0; c = (a+b)/2 S = S(a,c)+S(c,b) e = (S - S(a,b))/15 [V, TotalError] = AdaptQuad(a,b,S,e,V,TotalError) [V, TotalError] = function AdaptQuad(low,up,S,e,V,TotalError) if ( |e| >= tol*(up-low)/(b-a) ) mid = (low+up)/2 clow = (low+mid)/2 Slow = S(low,clow)+S(clow,mid) elow = (S-S(low,mid))/15 [V,TotalError] = AdaptQuad(low,mid,Slow,elow,V,TotalError] cup = (mid+up)/2 Sup = S(mid,cup)+S(cup,up) eup = (S-S(mid,up))/15 [V,TotalError] = AdaptQuad(mid,up,Sup,eup,V,TotalError] else V = V + S TotalError = TotalError + e endif It is possible to optimize this to avoid redundant evaluations of f(x). Matlab's quad function implements this (type "help quad" in Matlab). Matlab will plot the values at which f(x) is evaluated, as the integration proceeds. Ex: q = quad('sin',0,2*pi,tol,1) for tol = 1e-3, 1e-4, 1e-5 (true value = 0) q = quad('intfun',-1,1,tol,1) for tol = 1e-3, 1e-4, 1e-5 intfun(x) = min( 1/abs(x), exp(4) ), (true value = 10) q = quad('rsqrt',1e-20,1,tol,1) for tol = 1e-3, 1e-4, 1e-5 rsqrt(x) = 1/sqrt(x) (true value = 2-2e-10) This kind of adaptive approach can be used for approximation too, as discussed in section 6.14 of the text. Romberg Integration, aka Richardson Extrapolation Now we consider the case where we know the function is quite smooth, in which case we know there are high order derivatives, and we can expect a very accuracy solution by using a high order polynomial interpolant. When we studied differentiation, we saw how to take appropriate linear combinations of approximate derivatives D1(h) and D1(h/2), with step sizes h and h/2 and truncation errors O(h^2), to get a new approximation D2(h) with smaller truncation error O(h^4). We used a similar idea for adaptive quadrature. Since the trick worked well twice, let's use it again ("Any trick used at least 3 times is a standard technique") Let Tr{1}(h) be the value of the integral using the composite Trapezoidal rule with step size h. We saw that it has error -f''(z)*h^2/12 = O(h^2). Then we can let Tr{2}(h) = Tr{1}(h/2) + (Tr{1}(h/2) - Tr{1}(h))/4 and get a method with error O(h^4). Continuing recursively, we get (*) Tr{i+1}(h) = Tr{i}(h/2) + (Tr{i}(h/2) - Tr{i}(h))/(4^i -1) with error O(h^(2*i+2)). Note that each application of (*) requires us to evaluate the Trapezoidal rule at h/2, i.e. twice as many points as as h. But we have already evaluated f(.) at half these points in the course of evaluating the Trapezoidal rule at h. Still the number of points grows like 2^i, where i is the number of times we apply (*). The correctness of this depends on having an expansion of the form Tr{1}(h) = correct result of the integral + sum (from j=2 to infinity step 2) c{j}*h^(2*j) in place of the Taylor expansion, which we used for the derivative. (If the sum goes to a finite k instead of infinity, then we can only repeat Richardson Extrapolation k/2-1 times. The justification for this expansion is the Euler-Maclaurin Formula, which we will state but not prove (see sec 7.7 of book): integral (from 0 to 1) f(x)dx = Tr{1}(h) + sum (from k=1 to m-1) B(2*k)/(2*k)!*h^(2*k) * (f^(2*k-1)(1)-f^(2*k-1)(0)) - B(2*m)/(2*m)!*h^(2*m+1)*f^(2*m)(z) where B(2*m) is called a Bernoulli number (a constant). The idea of the proof is to repeatedly integrate by parts. Gaussian Quadrature: This is probably the best method for very smooth integrands, when only modest accuracy is required, because the number of points required can be anything, not just a power of 2 as with Romberg integration. In the Method of Undetermined Coefficients, we pointed out that using polynomial interpolation with a polynomial of degree n leads to a formula integral (from a to b) f(x) dx = sum {from i=0 to m} A{i}*f(x{i}) where the x{i} can be any n+1 distinct points between a and b, and the A{i} can be chosen to exactly integrate polynomials of degree n. The question arises as to what the best choice of x{i} is, and whether we can do better than degree n polynomials. The answers is yes: Theorem. x{0},...,x{n} and A{0},...,A{n} can be chosen to exactly integrate polynomials of degree 2*n+1. More generally, given any fixed positive weight function w(x), they can be chosen so that integral (from a to b) f(x)*w(x)dx = sum (from i=0 to n) A{i}*f(x{i}) exactly for all polynomials f(x) of degree at most 2*n+1. This choice of x{i} and A{i} is called Gaussian Quadrature. It is probably the fastest method for smooth integrands and modest accuracy requirements. Ex: integral (from -1 to 1) f(x)dx = sum (from i=0 to n) A{i}*f(x{i}) exactly for degree 2*n+1 polynomials if x{i} = cos((2*i+1)/(2*n+2)*pi) A{i} = pi/(n+1)*sin((2*i+1)/(2*n+2)*pi) In general, it is a lot of work to compute the x{i} and A{i} for given n and w(x), but they can be precomputed and stored in tables. Their computation is an exercise in linear algebra, to which we will return later. Finding integration software Within Matlab, we have already discussed quad(). For details, type "help quad". For programs in other programming languages, see gams.nist.gov, click on "Problem Decision Tree" then on "Differentiation, integration" then on "Quadrature" then further, depending on what one needs...