Numerical integration: Gaussian quadrature rules

1 APMA 0160 (A. Yew) Spring 2011 Numerical integration: Gaussian quadrature rulesMatlab s built-in Numerical integration function[Q,fcount]=quad(f,a,b,tol)is essentially oursimp_compextrcode with some further efficiency-enhancing thatquadrequires scalar functions to be defined with elementwise operations, sof(x) =21+x2should be entered asf=inline( (1+x.^2) , x )orf=@(x) (1+x.^2)The default tolerance forquadis 10 has another efficient integration command calledquadl, with the same input and outputarguments. The method underlyingquadlis a Gaussian quadrature rule.

2 Recall that each Newton Cotes quadrature rule came from integrating the Lagrange polynomial thatinterpolates the integrandfatnequally spacednodes in the interval [a,b]. Thus, in general, we expectthe degree of exactness of the rule to ben 1 (though, as we ve seen, some rules turn out to have ahigher-than-expected degree of exactness ).Is it possible to find quadrature rules that usennodes but have degree of exactness higher thann 1?Any (non-composite) quadrature rule can be written in the form(?)n i=1wif(xi)where thexiare nodes and thewiare non-zero constants called quadrature weights.

3 If thennodes are required to be equally spaced, then once you are told whether or not to include theendpoints ( whether the rule should be closed or open), the positions of the nodes are allfixed, soit only remains to find thenweightswi. These are given bywi= baLi(x)dxwhereLi(x) =n k=1, k6=ix xixk xiBut what if the nodes are not required to be equally spaced? If we put no restrictions on the positions ofthe nodes, other than that they must lie in the interval [a,b], then (?) would contain 2nfree parameters:w1,..,wnandx1.

4 ,xn. In general, 2nparameters could potentially define a polynomial of degreeup to 2n 1. Therefore, the highest degree of exactness we can expect to achieve with any quadraturerule is 2n illustrate the procedure of finding such a quadrature rule with degree of exactness 2n 1, let usconsider how to choose thewiandxiwhenn= 2 and the interval of integration is [ 1,1]. In otherwords, we want to determinew1,w2andx1,x2such thatw1f(x1) +w2f(x2) gives the exact value of 1 1f(x)dxwheneverfis a polynomial of degree 2(2) 1 = 3 or can be proved using the theory of orthogonal polynomials that, on the integration interval [ 1,1],to achieve the maximum degree of exactness 2n 1 withnnodes andnweights we should choose thenodesx1.

5 ,xnto be the roots (zeros) of the degree-nLegendre polynomialPn(x); the weights arethen given bywi= 1 1n k=1, k6=ix xixk xidx, and son i=1wif(xi) is an approximation of 1 1f(x) Legendre polynomials can be defined via the recursive relationPk+1(x) =2k+ 1k+ 1xPk(x) kk+ 1Pk 1(x) withP0(x) = 1,P1(x) = first few areP2(x) =32(x2 13),P3(x) =52(x3 35x) andP4(x) =18(35x4 30x2+ 3).EachPk(x) is a polynomial of degreek, and haskroots that all lie in the interval ( 1,1).So, to find the quadrature rule with maximum degree of exactness usingnnodes andnweights, inprinciple we need to: Find the Legendre polynomialPn(x).

6 Find the roots ofPn(x) in ( 1,1); these will be our nodesx1,..,xn(and so the quadraturerule will be an open rule). Find the Lagrange polynomial that interpolates the integrandf(x) atx1,..,xn. Integrate the Lagrange polynomial to determine the weightsw1,.., , the roots of the Legendre polynomials and their corresponding weights have been exten-sively tabulated, so we can simply use these tables without redoing the calculations. The only issueto be careful about is that the tabulated values were computed by taking [ 1,1] as the integrationinterval, whereas in any given problem the integration interval may not be [ 1,1].

7 Therefore we needto transform the tabulated values to analogous values on a general interval [a,b]. This is done throughthe following change of variable: 1 x 1 a x bif x ( 1)1 ( 1)=x ab aso thatx=b a2 x+a+b2 baf(x)dx= 1 1f(b a2 x+a+b2)dxd xd x= 1 1f(b a2 x+a+b2)b a2d x=b a2 1 1f(b a2 x+a+b2)d x b a2n i=1wif(b a2 xi+a+b2)If we defineh=b a(the length of the interval) andc=a+b2(the midpoint of the interval), thenthe roots xiin [ 1,1] are transformed to the nodesxiin [a,b] viaxi=h2 xi+c, and the quadratureformula for approximating baf(x)dxwill beh2timesthe formula for approximating the equivalentintegral over [ 1,1].

8 The quadrature rules defined above, using the roots of Legendre polynomials as their nodes, are calledGauss Legendre rules . They have degree of exactness 2n 1 (and order 2n). Gauss Legendrerules are open rules , and because the nodes are often positioned at irrational points in the interval,when we code the adaptive composite rules by repeatedly halving the interval, many extra functionevaluations may need to be order to make use of previously computed function values more efficiently, one could try to define aclosedGaussian quadrature rule.

9 Such a rule would havex1=aandxn=b, and it turns out that theappropriate choice of then 2 interior nodes should be the (transformed) roots ofP n 1(x) in ( 1,1).These roots and their associated weights are also available in tables, and the same transformation asabove would convert them to the corresponding nodes and quadrature formula on [a,b]. The rulesdefined in this way are calledGauss Lobatto rules . Note that becausex1andxnare now fixed,the formula ni=1wif(xi) has only 2n 2 free parameters, so the degree of exactness attainable by aGauss Lobatto rule is 2n 3 (order 2n 2).

10 Gauss Legendre nodes and coefficientsnRoots ofPn(x)Coefficientsci2 1/ 311/ 313 3/55/908/9 3/55/94 (15 + 2 30)/35(18 30)/36 (15 2 30)/35(18 + 30)/36 (15 2 30)/35(18 + 30)/36 (15 + 2 30)/35(18 30)/36 Gauss Lobatto nodes and coefficientsnx1, roots ofP n 1(x),xnCoefficientsci2 11113 11/304/311/34 11/6 1/ 55/61/ 55/611/65 11/10 3/749/90032/45 3/749/9011/10 Matlab squadlemploys a composite 4-node Gauss Lobatto rule with some further efficiency-enhancingfeatures.