Chapter 6 Quadrature - MathWorks

1 Chapter 6 QuadratureThe termnumerical integrationcovers several different tasks, including numericalevaluation of integrals and numerical solution of ordinary differential equations. Sowe use the somewhat old-fashioned termquadraturefor the simplest of these, thenumerical evaluation of a definite integral. Modern Quadrature algorithms auto-matically vary an adaptive step Adaptive QuadratureLetf(x) be a real-valued function of a real variable, defined on a finite intervala x b. We seek to compute the value of the integral, baf(x)dx:The word Quadrature reminds us of an elementary technique for finding thisarea plot the function on graph paper and count the number of little squares thatlie underneath the Figure , there are 148 little squares underneath the curve.

2 If the areaof one little square is 3=512, then a rough estimate of the integral is 148 3=512 =0 quadratureinvolves careful selection of the points wheref(x) is sam-pled. We want to evaluate the function at as few points as possible while approx-imating the integral to within some specified accuracy. A fundamental additiveproperty of a definite integral is the basis for adaptive Quadrature . Ifcis any pointbetweenaandb, then baf(x)dx= caf(x)dx+ bcf(x)dx:The idea is that if we can approximate each of the two integrals on the right towithin a specified tolerance, then the sum gives us the desired result. If not, weSeptember 16, 201312 Chapter 6. QuadratureFigure recursively apply the additive property to each of the intervals [a;c] and [c;b].The resulting algorithm will adapt to the integrand automatically, partitioning theinterval into subintervals with fine spacing where the integrand is varying rapidlyand coarse spacing where the integrand is varying Basic Quadrature RulesThe derivation of the Quadrature rule used by ourMatlabfunction begins withtwo of the basic Quadrature rules shown in Figure : themidpoint ruleand thetrapezoid rule.

3 Leth=b abe the length of the interval. The midpoint rule,M,approximates the integral by the area of a rectangle whose base has lengthhandwhose height is the value of the integrand at the midpoint:M=hf(a+b2):The trapezoid rule,T, approximates the integral by the area of a trapezoid withbasehand sides equal to the values of the integrand at the two endpoints:T=hf(a) +f(b)2:The accuracy of a Quadrature rule can be predicted in part by examining itsbehavior on polynomials. Theorderof a Quadrature rule is the degree of the lowestdegree polynomial that the rule does not integrate exactly. If a Quadrature rule oforderpis used to integrate a smooth function over a small interval of lengthh, thena Taylor series analysis shows that the error is proportional tohp. The Basic Quadrature Rules3 Midpoint ruleTrapezoid ruleSimpson s ruleComposite Simpson s ruleFigure Quadrature and the trapezoid rule are both exact for constant and linear functions ofx,but neither of them is exact for a quadratic inx, so they both have order two.

4 (Theorder of a rectangle rule with heightf(a) orf(b) instead of the midpoint is onlyone.)The accuracy of the two rules can be compared by examining their behavioron the simple integral 10x2dx=13:The midpoint rule givesM= 1(12)2=14:The trapezoid rule givesT= 1(0 + 122)=12:So the error inMis 1=12, while the error inTis 1=6. The errors have oppo-site signs and, perhaps surprisingly, the midpoint rule is twice as accurate as thetrapezoid 6. QuadratureThis turns out to be true more generally. For integrating smooth functionsover short intervals,Mis roughly twice as accurate asTand the errors have oppositesigns. Knowing these error estimates allows us to combine the two and get a rulethat is usually more accurate than either one separately. If the error inTwereexactly 2 times the error inM, then solvingS T= 2(S M)forSwould give us the exact value of the integral.

5 In any case, the solutionS=23M+13 Tis usually a more accurate approximation than eitherMorTalone. This ruleis known asSimpson's rule. It can also be derived by integrating the quadraticfunction that interpolates the integrand at the two endpoints,aandb, and themidpoint,c= (a+b)=2:S=h6(f(a) + 4f(c) +f(b)):It turns out thatSalso integrates cubics exactly, but not quartics, so its orderis can carry this process one step further using the two halves of the interval,[a;c] and [c;b]. Letdandebe the midpoints of these two subintervals:d= (a+c)=2ande= (c+b)=2. Apply Simpson s rule to each subinterval to obtain a quadraturerule over [a;b]:S2=h12(f(a) + 4f(d) + 2f(c) + 4f(e) +f(b)):This is an example of acompositequadrature rule. See Figure the same integral, so their difference can be used as anestimate of the error:E= (S2 S):Moreover, the two can be combined to get an even more accurate approximation,Q.

6 Both rules are of order four, but theS2step size is half theSstep size, soS2isroughly 24times as accurate. Thus,Qis obtained by solvingQ S= 16(Q S2):The result isQ=S2+ (S2 S)=15:Exercise asks you to expressQas a weighted combination of the five functionvaluesf(a) throughf(e) and to establish that its order is six. The rule is knownasWeddle's rule, the sixth-orderNewton{Cotes rule, and also as the first step ofRomberg integration . We will simply call it theextrapolated Simpson's rulebecauseit uses Simpson s rule for two different values ofhand then extrapolates towardh= quadtx, quadtx, quadguiTheMatlabfunctionquaduses the extrapolated Simpson s rule in an adaptiverecursive algorithm. Our textbook functionquadtxis a simplified version functionquadguiprovides a graphical demonstration of the behavior ofquadandquadtx.}

7 It produces a dynamic plot of the function values selected by theadaptive algorithm. The count of function evaluations is shown in the title positionon the initial portion ofquadtxevaluates the integrandf(x) three times togive the first, unextrapolated, Simpson s rule estimate. A recursive subfunction,quadtxstep, is then called to complete the [Q,fcount] = quadtx(F,a,b,tol,varargin)%QUADTX Evaluate definite integral Q = QUADTX(F,A,B) approximates the integral of F(x)% from A to B to within a tolerance of Q = QUADTX(F,A,B,tol) uses tol instead of The first argument, F, is a function handle or% an anonymous function that defines F(x).%% Arguments beyond the first four,% Q = QUADTX(F,a,b,tol,p1,p2,..), are passed on to the% integrand, F(x,p1,p2.)

8 %% [Q,fcount] = QUADTX(F,..) also counts the number of% evaluations of F(x).%% See also QUAD, QUADL, DBLQUAD, Default toleranceif nargin < 4 | isempty(tol)tol = ;end% Initializationc = (a + b)/2;fa = F(a,varargin{:});fc = F(c,varargin{:});fb = F(b,varargin{:});% Recursive call[Q,k] = quadtxstep(F, a, b, tol, fa, fc, fb, varargin{:});fcount = k + 3;6 Chapter 6. QuadratureEach recursive call ofquadtxstepcombines three previously computed func-tion values with two more to obtain the two Simpson s approximations for a par-ticular interval. If their difference is small enough, they are combined to return theextrapolated approximation for that interval. If their difference is larger than thetolerance, the recursion proceeds on each of the two half [Q,fcount] = quadtxstep(F,a,b,tol,fa,fc,fb,varargin)% Recursive subfunction used by = b - a;c = (a + b)/2;fd = F((a+c)/2,varargin{:});fe = F((c+b)/2,varargin{:});Q1 = h/6 * (fa + 4*fc + fb);Q2 = h/12 * (fa + 4*fd + 2*fc + 4*fe + fb);if abs(Q2 - Q1) <= tolQ = Q2 + (Q2 - Q1)/15;fcount = 2;else[Qa,ka] = quadtxstep(F, a, c, tol, fa, fd, fc, varargin{:});[Qb,kb] = quadtxstep(F, c, b, tol, fc, fe, fb, varargin{:});Q = Qa + Qb;fcount = ka + kb + 2;endThe choice of tolerance for comparison with the error estimates is important,but a little tricky.

9 If a tolerance is not specified as the fourth argument to thefunction, then 10 6is used as the tricky part is how to specify the tolerance in the recursive calls. Howsmall does the tolerance in each recursive call have to be in order for the final resultto have the desired accuracy? One approach would cut the tolerance in half witheach level in the recursion. The idea is that if bothQaandQbhave errors less thantol/2, then their sum certainly has an error less thantol. If we did this, the twostatements[Qa,ka] = quadtxstep(F, a, c, tol, fa, fd, fc, varargin{:});[Qb,kb] = quadtxstep(F, c, b, tol, fc, fe, fb, varargin{:});would havetol/2in place , this approach is too conservative. We are estimating the error inthe two separate Simpson s rules, not their extrapolated combination.

10 So the actualerror is almost always much less than the estimate. More importantly, the subinter-vals where the actual error is close to the estimate are usually fairly rare. We canallow one of the two recursive calls to have an error close to the tolerance becausethe other subinterval will probably have a much smaller error. For these reasons,the same value oftolis used in each recursive Specifying Integrands7 Our textbook function does have one serious defect: there is no provision forfailure. It is possible to try to evaluate integrals that do not exist. For example, 1013x 1dxhas a nonintegrable singularity. Attempting to evaluate this integral withquadtxresults in a computation that runs for a long time and eventually terminates withan error message about the maximum recursion limit.