Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 40

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 40)

. . , xN , you need to find the weights wj ,j = 1, . . . , N . One way to do this (not the most efficient) is to solve the set oflinear equations p (x )0 1 p1 (x1 )...pN−1 (x1 )...... 'bw1a W (x)p0 (x)dx  w2  0 .  = ..  ..  .p0 (xN )p1 (xN ).... . . pN−1 (xN )wN(4.5.8)0Equation (4.5.8) simply solves for those weights such that the quadrature (4.5.1)gives the correct answer for the integral of the first N orthogonal polynomials. Notethat the zeros on the right-hand side of (4.5.8) appear because p1 (x), . .

. , pN−1 (x)are all orthogonal to p0 (x), which is a constant. It can be shown that, with thoseweights, the integral of the next N − 1 polynomials is also exact, so that thequadrature is exact for all polynomials of degree 2N − 1 or less. Another way toevaluate the weights (though one whose proof is beyond our scope) is by the formulawj =pN−1 |pN−1 pN−1 (xj )pN (xj )(4.5.9)where pN (xj ) is the derivative of the orthogonal polynomial at its zero xj .The computation of Gaussian quadrature rules thus involves two distinct phases:(i) the generation of the orthogonal polynomials p0 , .

. . , pN , i.e., the computation ofthe coefficients aj , bj in (4.5.6); (ii) the determination of the zeros of pN (x), andthe computation of the associated weights. For the case of the “classical” orthogonalpolynomials, the coefficients aj and bj are explicitly known (equations 4.5.10 –4.5.14 below) and phase (i) can be omitted. However, if you are confronted with a“nonclassical” weight function W (x), and you don’t know the coefficients aj andbj , the construction of the associated set of orthogonal polynomials is not trivial.We discuss it at the end of this section.Computation of the Abscissas and WeightsThis task can range from easy to difficult, depending on how much you alreadyknow about your weight function and its associated polynomials.

In the case ofclassical, well-studied, orthogonal polynomials, practically everything is known,including good approximations for their zeros. These can be used as starting guesses,enabling Newton’s method (to be discussed in §9.4) to converge very rapidly.Newton’s method requires the derivative pN (x), which is evaluated by standardrelations in terms of pN and pN−1 .

The weights are then conveniently evaluated byequation (4.5.9). For the following named cases, this direct root-finding is faster,by a factor of 3 to 5, than any other method.Here are the weight functions, intervals, and recurrence relations that generatethe most commonly used orthogonal polynomials and their corresponding Gaussianquadrature formulas.4.5 Gaussian Quadratures and Orthogonal Polynomials151Gauss-Legendre:−1<x<1W (x) = 1(j + 1)Pj+1 = (2j + 1)xPj − jPj−1(4.5.10)Gauss-Chebyshev:W (x) = (1 − x2 )−1/2−1<x<1Tj+1 = 2xTj − Tj−1(4.5.11)Gauss-Laguerre:W (x) = xα e−x0<x<∞αα(j + 1)Lαj+1 = (−x + 2j + α + 1)Lj − (j + α)Lj−1(4.5.12)Gauss-Hermite:W (x) = e−x2−∞<x<∞Hj+1 = 2xHj − 2jHj−1(4.5.13)Gauss-Jacobi:W (x) = (1 − x)α (1 + x)β(α,β)(α,β)cj Pj+1 = (dj + ej x)Pj−1<x<1(α,β)− fj Pj−1(4.5.14)where the coefficients cj , dj , ej , and fj are given bycj = 2(j + 1)(j + α + β + 1)(2j + α + β)dj = (2j + α + β + 1)(α2 − β 2 )ej = (2j + α + β)(2j + α + β + 1)(2j + α + β + 2)(4.5.15)fj = 2(j + α)(j + β)(2j + α + β + 2)We now give individual routines that calculate the abscissas and weights forthese cases.

First comes the most common set of abscissas and weights, those ofGauss-Legendre. The routine, due to G.B. Rybicki, uses equation (4.5.9) in thespecial form for the Gauss-Legendre case,wj =2(1 −x2j )[PN (xj )]2(4.5.16)The routine also scales the range of integration from (x1 , x2) to (−1, 1), and providesabscissas xj and weights wj for the Gaussian formula& x2Nf(x)dx =wj f(xj )(4.5.17)x1j=1152Chapter 4.#include <math.h>#define EPS 3.0e-11Integration of FunctionsEPS is the relative precision.void gauleg(float x1, float x2, float x[], float w[], int n)Given the lower and upper limits of integration x1 and x2, and given n, this routine returnsarrays x[1..n] and w[1..n] of length n, containing the abscissas and weights of the GaussLegendre n-point quadrature formula.{int m,j,i;double z1,z,xm,xl,pp,p3,p2,p1;High precision is a good idea for this routine.m=(n+1)/2;The roots are symmetric in the interval, soxm=0.5*(x2+x1);we only have to find half of them.xl=0.5*(x2-x1);for (i=1;i<=m;i++) {Loop over the desired roots.z=cos(3.141592654*(i-0.25)/(n+0.5));Starting with the above approximation to the ith root, we enter the main loop ofrefinement by Newton’s method.do {p1=1.0;p2=0.0;for (j=1;j<=n;j++) {Loop up the recurrence relation to get thep3=p2;Legendre polynomial evaluated at z.p2=p1;p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;}p1 is now the desired Legendre polynomial.

We next compute pp, its derivative,by a standard relation involving also p2, the polynomial of one lower order.pp=n*(z*p1-p2)/(z*z-1.0);z1=z;z=z1-p1/pp;Newton’s method.} while (fabs(z-z1) > EPS);x[i]=xm-xl*z;Scale the root to the desired interval,x[n+1-i]=xm+xl*z;and put in its symmetric counterpart.w[i]=2.0*xl/((1.0-z*z)*pp*pp);Compute the weightw[n+1-i]=w[i];and its symmetric counterpart.}}Next we give three routines that use initial approximations for the roots givenby Stroud and Secrest [2].

The first is for Gauss-Laguerre abscissas and weights, tobe used with the integration formula&0#include <math.h>#define EPS 3.0e-14#define MAXIT 10∞xα e−x f(x)dx =Nwj f(xj )(4.5.18)j=1Increase EPS if you don’t have this precision.void gaulag(float x[], float w[], int n, float alf)Given alf, the parameter α of the Laguerre polynomials, this routine returns arrays x[1..n]and w[1..n] containing the abscissas and weights of the n-point Gauss-Laguerre quadratureformula.

The smallest abscissa is returned in x[1], the largest in x[n].{float gammln(float xx);void nrerror(char error_text[]);int i,its,j;float ai;4.5 Gaussian Quadratures and Orthogonal Polynomials153High precision is a good idea for this routine.for (i=1;i<=n;i++) {Loop over the desired roots.if (i == 1) {Initial guess for the smallest root.z=(1.0+alf)*(3.0+0.92*alf)/(1.0+2.4*n+1.8*alf);} else if (i == 2) {Initial guess for the second root.z += (15.0+6.25*alf)/(1.0+0.9*alf+2.5*n);} else {Initial guess for the other roots.ai=i-2;z += ((1.0+2.55*ai)/(1.9*ai)+1.26*ai*alf/(1.0+3.5*ai))*(z-x[i-2])/(1.0+0.3*alf);}for (its=1;its<=MAXIT;its++) {Refinement by Newton’s method.p1=1.0;p2=0.0;for (j=1;j<=n;j++) {Loop up the recurrence relation to get thep3=p2;Laguerre polynomial evaluated at z.p2=p1;p1=((2*j-1+alf-z)*p2-(j-1+alf)*p3)/j;}p1 is now the desired Laguerre polynomial.

We next compute pp, its derivative,by a standard relation involving also p2, the polynomial of one lower order.pp=(n*p1-(n+alf)*p2)/z;z1=z;z=z1-p1/pp;Newton’s formula.if (fabs(z-z1) <= EPS) break;}if (its > MAXIT) nrerror("too many iterations in gaulag");x[i]=z;Store the root and the weight.w[i] = -exp(gammln(alf+n)-gammln((float)n))/(pp*n*p2);}double p1,p2,p3,pp,z,z1;}Next is a routine for Gauss-Hermite abscissas and weights. If we use the“standard” normalization of these functions, as given in equation (4.5.13), we findthat the computations overflow for large N because of various factorials that occur.We can avoid this by using instead the orthonormal set of polynomials Hj .

Theyare generated by the recurrence*)12jHj −Hj−1(4.5.19)H−1 = 0, H0 = 1/4 , Hj+1 = xj+1j +1πThe formula for the weights becomeswj =2(4.5.20)(Hj )2while the formula for the derivative with this normalization is(Hj = 2j Hj−1(4.5.21)The abscissas and weights returned by gauher are used with the integration formula&∞−∞e−x f(x)dx =2Nj=1wj f(xj )(4.5.22)154Chapter 4.#include <math.h>#define EPS 3.0e-14#define PIM4 0.7511255444649425#define MAXIT 10Integration of FunctionsRelative precision.1/π1/4 .Maximum iterations.void gauher(float x[], float w[], int n)Given n, this routine returns arrays x[1..n] and w[1..n] containing the abscissas and weightsof the n-point Gauss-Hermite quadrature formula. The largest abscissa is returned in x[1], themost negative in x[n].{void nrerror(char error_text[]);int i,its,j,m;double p1,p2,p3,pp,z,z1;High precision is a good idea for this routine.m=(n+1)/2;The roots are symmetric about the origin, so we have to find only half of them.for (i=1;i<=m;i++) {Loop over the desired roots.if (i == 1) {Initial guess for the largest root.z=sqrt((double)(2*n+1))-1.85575*pow((double)(2*n+1),-0.16667);} else if (i == 2) {Initial guess for the second largest root.z -= 1.14*pow((double)n,0.426)/z;} else if (i == 3) {Initial guess for the third largest root.z=1.86*z-0.86*x[1];} else if (i == 4) {Initial guess for the fourth largest root.z=1.91*z-0.91*x[2];} else {Initial guess for the other roots.z=2.0*z-x[i-2];}for (its=1;its<=MAXIT;its++) {Refinement by Newton’s method.p1=PIM4;p2=0.0;for (j=1;j<=n;j++) {Loop up the recurrence relation to getp3=p2;the Hermite polynomial evaluated atp2=p1;z.p1=z*sqrt(2.0/j)*p2-sqrt(((double)(j-1))/j)*p3;}p1 is now the desired Hermite polynomial.

Характеристики

Список файлов книги