6 Numerical Integration - Department of Mathematics

1 D. Levy6 Numerical Basic ConceptsIn this chapter we are going to explore various ways for approximating the integral of afunction over a given domain. There are various reasons as of why such approximationscan be useful. First, not every function can be analytically integrated. Second, even if aclosed Integration formula exists, it might still not be the most efficient way of calculatingthe integral. In addition, it can happen that we need to integrate an unknown function,in which only some samples of the function are order to gain some insight on Numerical Integration , it is natural to review Rie-mann Integration , a framework that can be viewed as an approach for approximat-ing integrals. We assume thatf(x) is a bounded function defined on [a, b] and that{x0, .. , xn}is a partition (P) of [a, b]. For eachiwe letMi(f) = supx [xi 1,xi]f(x),andmi(f) = infx [xi 1,xi]f(x),Letting xi=xi xi 1,the upper (Darboux) sumoff(x) with respect to thepartitionPis defined asU(f, P) =n i=1Mi xi,( )whilethe lower (Darboux) sumoff(x) with respect to the partitionPis defined asL(f, P) =n i=1mi xi.

2 ( )The upper integraloff(x) on [a, b] is defined asU(f) = inf(U(f, P)),andthe lower integraloff(x) is defined asL(f) = sup(L(f, P)),where both the infimum and the supremum are taken over all possible partitions,P, ofthe interval [a, b]. If the upper and lower integral off(x) are equal to each other, theircommon value is denoted by baf(x)dxand is referred to asthe Riemann integraloff(x). Basic ConceptsD. LevyFor the purpose of the present discussion we can think of the upper and lower Dar-boux sums ( ), ( ), as two approximations of the integral (assuming that the functionis indeed integrable). Of course, these sums are not defined in the most convenient wayfor an approximation algorithm. This is because we need to find the extrema of the func-tion in every subinterval. Finding the extrema of the function, may be a complicatedtask on its own, which we would like to simpler approach for approximating the value of baf(x)dxwould be to computethe product of the value of the function at one of the end-points of the interval by thelength of the interval.

3 In case we choose the end-point where the function is evaluatedto bex=a, we obtain baf(x)dx f(a)(b a).( )This approximation ( ) is calledthe rectangular method(see Figure ). Numeri-cal Integration formulas are also referred to asintegration rulesorquadratures, andhence we can refer to ( ) as the rectangular rule or the rectangular quadrature. Thepointsx0, .. xnthat are used in the quadrature formula are calledquadrature (a)f(b)xf(x)Figure : A rectangular quadratureA variation on the rectangular rule isthe midpoint rule. Similarly to the rectan-gular rule, we approximate the value of the integral baf(x)dxby multiplying the lengthof the interval by the value of the function at one point. Only this time, we replace thevalue of the function at an endpoint, by the value of the function at the center point2D. Basic Concepts12(a+b), , baf(x)dx (b a)f(a+b2).( )(see also Fig ). As we shall see below, the midpoint quadrature ( ) is a moreaccurate quadrature than the rectangular rule ( ).

4 A(a+b)/2bf(a)f((a+b)/2)f(b)xf(x)Figure : A midpoint quadratureIn order to compute the quadrature error for the midpoint rule ( ), we considerthe primitive functionF(x),F(x) = xaf(s)ds,and expand a+haf(s)ds=F(a+h) =F(a) +hF (a) +h22F (a) +h36F (a) +O(h4) ( )=hf(a) +h22f (a) +h36f (a) +O(h4)If we letb=a+h, we have (expandingf(a+h/2)) for the quadrature error,E,E= a+haf(s)ds hf(a+h2)=hf(a) +h22f (a) +h36f (a) +O(h4) h[f(a) +h2f (a) +h28f (a) +O(h3)], Basic ConceptsD. Levywhich means that the error term is of orderO(h3). Having an error of orderh3doesnot mean that this is a third-order method. In our case, the parameterhequals tob a. It is not a parameter that should be vied as a small value that goes to zero. It isfixed. The error of the midpoint method is of the order ofO((b a)3). Unfortunately,these calculations cannot directly provide us with an accurate estimate of the is the case since when truncating two Taylor approximations, we are left with anerror terms that are evaluated at two (generally different) intermediate points.

5 Hencewe cannot directly combine the error term16h3f ( 1) with 18h3f ( 2). This can still bedone, but we have to use a better main difficulty in evaluating the difference between the exact value, baf(x)dx,and its midpoint rule approximation, (b a)f(a+b2), is due to having an integral in oneterm and no integral in the second term. The approach will be to replace the midpointapproximation with an integral expression. Indeed, if we denote the midpoint byc, ,c=a+b2,then the tangent line tof(x) atx=cis given byP1(x) =f(c) +f (c)(x c).Clearly, baP1(x)dx= (b a)f(c),and hence baf(x)dx (b a)f(a+b2)= ba(f(x) P1(x)) estimate the difference betweenf(x) andP1(x) we can expandf(x) aroundx= thatx [a, b], we havef(x) =f(c+ (x c)) =f(c) +f (c)(x c) +12f ( )(x c)2, (a, b).Hence ba(f(x) P1(x))dx= ba12f ( x)(x c) view of the midvalue theorem for integrals, the last integral can be replaces by12f ( ) ba(x c)2dx=124(b a)3f ( ),a < < b.

6 ( ) this section we assumed that all functions we are interested inintegrating are actually integrable in the domain of interest. We also assumed that theyare bounded and that they are defined at every point, so that whenever we need toevaluate a function at a point, we can do it. We will go on and use these assumptionsthroughout the Integration via Integration via InterpolationOne direct way of obtaining quadratures from given samples of a function is by integrat-ing an interpolant. As always, our goal is to evaluateI= baf(x)dx. We assume thatthe values of the functionf(x) are given atn+ 1 points:x0, .. , xn [a, b]. Note thatwe do not require the first pointx0to be equal toa, and the same holds for the rightside of the interval. Given the valuesf(x0), .. f(xn), we can write the interpolatingpolynomial of degree6n, which in the Largenge form isPn(x) =n i=0f(xi)li(x),withli(x) =n j=0j6=ix xjxi xj, integral off(x) can then be approximated by the integral ofPn(x), , baf(x)dx baPn(x)dx=n i=0f(xi) bali(x)dx=n i=0 Aif(xi).

7 ( )The quadrature coefficientsAiin ( ) are given byAi= bali(x)dx.( )Note that if we want to integrate several different functions, and use their values atthe same points (x0, .. , xn), the quadrature coefficients ( ) should be computed onlyonce, since they do not depend on the function that is being integrated. If we change theinterpolation/ Integration points, then we must recompute the quadrature equally spaced points,x0, .. , xn, a Numerical Integration formula of the form baf(x)dx n i=0 Aif(xi),( )is called aNewton-Cotes letn= 1 and consider two interpolation points which we set asx0=a, x1= this casel0(x) =b xb a,l1(x) =x ab Integration via InterpolationD. LevyHenceA0= bal0(x) = bab xb adx=b ,A1= bal1(x) = bax ab adx=b a2= resulting quadrature is the so-calledtrapezoidal rule, badx b a2[f(a) +f(b)],( )(see Figure ).abf(a)f(b)xf(x)Figure : A trapezoidal quadratureWe can now use the interpolation error to compute the error in the quadrature ( ).

8 The interpolation error isf(x) P1(x) =12f ( x)(x a)(x b), x (a, b).We recall that according to the midvalue theorem for integrals, ifu(x) andv(x) arecontinuous on [a, b] and ifv>0, then there exists (a, b) such that bau(x)v(x)dx=u( ) bav(x) Integration via InterpolationHence, the interpolation error is given byE= ba12f ( x)(x a)(x b) =f ( )2 ba(x a)(x b)dx= f ( )12(b a)3,( )with (a, b). We note that the quadratures ( ),( ), are exact for polynomials of degree6n. For ifp(x) is a polynomial of degree6n, it can be written asp(x) =n i=0p(xi)li(x).(Two polynomials of degree6nthat agree with each other atn+ 1 points mustbe identical). Hence bap(x)dx=n i=0p(xi) bali(x)dx=n i=0 Aip(xi).2. As of the opposite direction. Assume that the quadrature baf(x)dx n i=0 Aif(xi),is exact for all polynomials of degree6n. We know thatdeg(lj(x)) =n,and hence balj(x)dx=n i=0 Ailj(xi) =n i=0Ai ij= menas that the quadrature coefficients must be given byAj= balj(x) Composite Integration RulesD.

9 Composite Integration RulesIn a composite quadrature, we divide the interval into subintervals and apply an inte-gration rule to each subinterval. We demonstrate this idea with a couple of the pointsa=x0< x1< < xn= trapezoidal ruleis obtained by applying the trapezoidal rule in eachsubinterval [xi 1, xi],i= 1, .. , n, , baf(x)dx=n i=1 xixi 1f(x)dx 12n i=1(xi xi 1)[f(xi 1) +f(xi)],( )(see Figure ).x0x1x2xn 1xnxf(x)Figure : A composite trapezoidal ruleA particular case is when these points are uniformly spaced, , when all intervalshave an equal length. For example, ifxi=a+ih,whereh=b an,8D. Composite Integration Rulesthen baf(x)dx h2[f(a) + 2n 1 i=1f(a+ih) +f(b)]=hn i=0 f(a+ih).( )The notation of a sum with two primes, , means that we sum over all the terms withthe exception of the first and last terms that are being divided by can also compute the error term as a function of the distance between neighboringpoints,h.

10 We know from ( ) that in every subinterval the quadrature error is h312f ( x).Hence, the overall error is obtained by summing overnsuch terms:n i=1 h312f ( i) = h3n12[1nn i=1f ( i)].Here, we use the notation ito denote an intermediate point that belongs to theithinterval. LetM=1nn i=1f ( i).Clearlyminx [a,b]f (x)6M6maxx [a,b]f (x)If we assume thatf (x) is continuous in [a, b] (which we anyhow do in order for theinterpolation error formula to be valid) then there exists a point [a, b] such thatf ( ) = (recalling that (b a)/n=h, we haveE= (b a)h212f ( ), [a, b].( )This means that the composite trapezoidal rule is second-order the interval [a, b] we assumensubintervals and leth=b Additional Integration TechniquesD. LevyThe quadrature points arexj=a+(j 12)h,j= 1,2, .. , composite midpoint rule is given by applying the midpoint rule ( ) in everysubinterval, , baf(x)dx hn j=1f(xj).( )Equation ( ) is known asthe composite midpoint order to obtain the quadrature error in the approximation ( ) we recall thatin each subinterval the error is given according to ( ), ,Ej=h324f ( j), j (xj h2, xj+h2).)