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

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

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

We next compute pp, its derivative, bythe relation (4.5.21) using p2, the polynomial of one lower order.pp=sqrt((double)2*n)*p2;z1=z;z=z1-p1/pp;Newton’s formula.if (fabs(z-z1) <= EPS) break;}if (its > MAXIT) nrerror("too many iterations in gauher");x[i]=z;Store the rootx[n+1-i] = -z;and its symmetric counterpart.w[i]=2.0/(pp*pp);Compute the weightw[n+1-i]=w[i];and its symmetric counterpart.}}Finally, here is a routine for Gauss-Jacobi abscissas and weights, whichimplement the integration formula&1−1(1 − x)α (1 + x)β f(x)dx =Nj=1wj f(xj )(4.5.23)4.5 Gaussian Quadratures and Orthogonal Polynomials#include <math.h>#define EPS 3.0e-14#define MAXIT 10155Increase EPS if you don’t have this precision.void gaujac(float x[], float w[], int n, float alf, float bet)Given alf and bet, the parameters α and β of the Jacobi polynomials, this routine returnsarrays x[1..n] and w[1..n] containing the abscissas and weights of the n-point Gauss-Jacobiquadrature formula.

The largest abscissa is returned in x[1], the smallest in x[n].{float gammln(float xx);void nrerror(char error_text[]);int i,its,j;float alfbet,an,bn,r1,r2,r3;double a,b,c,p1,p2,p3,pp,temp,z,z1;High 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 largest root.an=alf/n;bn=bet/n;r1=(1.0+alf)*(2.78/(4.0+n*n)+0.768*an/n);r2=1.0+1.48*an+0.96*bn+0.452*an*an+0.83*an*bn;z=1.0-r1/r2;} else if (i == 2) {Initial guess for the second largest root.r1=(4.1+alf)/((1.0+alf)*(1.0+0.156*alf));r2=1.0+0.06*(n-8.0)*(1.0+0.12*alf)/n;r3=1.0+0.012*bet*(1.0+0.25*fabs(alf))/n;z -= (1.0-z)*r1*r2*r3;} else if (i == 3) {Initial guess for the third largest root.r1=(1.67+0.28*alf)/(1.0+0.37*alf);r2=1.0+0.22*(n-8.0)/n;r3=1.0+8.0*bet/((6.28+bet)*n*n);z -= (x[1]-z)*r1*r2*r3;} else if (i == n-1) {Initial guess for the second smallest root.r1=(1.0+0.235*bet)/(0.766+0.119*bet);r2=1.0/(1.0+0.639*(n-4.0)/(1.0+0.71*(n-4.0)));r3=1.0/(1.0+20.0*alf/((7.5+alf)*n*n));z += (z-x[n-3])*r1*r2*r3;} else if (i == n) {Initial guess for the smallest root.r1=(1.0+0.37*bet)/(1.67+0.28*bet);r2=1.0/(1.0+0.22*(n-8.0)/n);r3=1.0/(1.0+8.0*alf/((6.28+alf)*n*n));z += (z-x[n-2])*r1*r2*r3;} else {Initial guess for the other roots.z=3.0*x[i-1]-3.0*x[i-2]+x[i-3];}alfbet=alf+bet;for (its=1;its<=MAXIT;its++) {Refinement by Newton’s method.temp=2.0+alfbet;Start the recurrence with P0 and P1 to avoidp1=(alf-bet+temp*z)/2.0;a division by zero when α + β = 0 orp2=1.0;−1.for (j=2;j<=n;j++) {Loop up the recurrence relation to get thep3=p2;Jacobi polynomial evaluated at z.p2=p1;temp=2*j+alfbet;a=2*j*(j+alfbet)*(temp-2.0);b=(temp-1.0)*(alf*alf-bet*bet+temp*(temp-2.0)*z);c=2.0*(j-1+alf)*(j-1+bet)*temp;p1=(b*p2-c*p3)/a;}pp=(n*(alf-bet-temp*z)*p1+2.0*(n+alf)*(n+bet)*p2)/(temp*(1.0-z*z));p1 is now the desired Jacobi polynomial.

We next compute pp, its derivative, bya standard relation involving also p2, the polynomial of one lower order.z1=z;z=z1-p1/pp;Newton’s formula.156Chapter 4.Integration of Functionsif (fabs(z-z1) <= EPS) break;}if (its > MAXIT) nrerror("too many iterations in gaujac");x[i]=z;Store the root and the weight.w[i]=exp(gammln(alf+n)+gammln(bet+n)-gammln(n+1.0)gammln(n+alfbet+1.0))*temp*pow(2.0,alfbet)/(pp*p2);}}Legendre polynomials are special cases of Jacobi polynomials with α = β = 0,but it is worth having the separate routine for them, gauleg, given above.

Chebyshevpolynomials correspond to α = β = −1/2 (see §5.8). They have analytic abscissasand weights:π(j − 12 )xj = cosN(4.5.24)πwj =NCase of Known RecurrencesTurn now to the case where you do not know good initial guesses for the zeros of yourorthogonal polynomials, but you do have available the coefficients aj and bj that generatethem. As we have seen, the zeros of pN (x) are the abscissas for the N -point Gaussianquadrature formula. The most useful computational formula for the weights is equation(4.5.9) above, since the derivative pN can be efficiently computed by the derivative of (4.5.6)in the general case, or by special relations for the classical polynomials. Note that (4.5.9) isvalid as written only for monic polynomials; for other normalizations, there is an extra factorof λN /λN −1 , where λN is the coefficient of xN in pN .Except in those special cases already discussed, the best way to find the abscissas is notto use a root-finding method like Newton’s method on pN (x).

Rather, it is generally fasterto use the Golub-Welsch [3] algorithm, which is based on a result of Wilf [4]. This algorithmnotes that if you bring the term xpj to the left-hand side of (4.5.6) and the term pj+1 to theright-hand side, the recurrence relation can be written in matrix form as   p0p0a0 10 p1   b1 a1 1  p1   0  .    .   . .... · . + . x ..  =   .   . ..    0 bN −2 aN −21pN −2pN −2pNbN −1 aN −1pN −1pN −1orxp = T · p + pN eN −1(4.5.25)Here T is a tridiagonal matrix, p is a column vector of p0 , p1 , .

. . , pN −1 , and eN −1 is a unitvector with a 1 in the (N − 1)st (last) position and zeros elsewhere. The matrix T can besymmetrized by a diagonal similarity transformation D to give√b1√a0√b2 b1 a1−1..J = DTD = (4.5.26)....√√bN −2 √aN −2bN −1bN −1 aN −1The matrix J is called the Jacobi matrix (not to be confused with other matrices namedafter Jacobi that arise in completely different problems!).

Now we see from (4.5.25) that1574.5 Gaussian Quadratures and Orthogonal PolynomialspN (xj ) = 0 is equivalent to xj being an eigenvalue of T. Since eigenvalues are preservedby a similarity transformation, xj is an eigenvalue of the symmetric tridiagonal matrix J.Moreover, Wilf [4] shows that if vj is the eigenvector corresponding to the eigenvalue xj ,normalized so that v · v = 1, then2wj = µ0 vj,1where&(4.5.27)bµ0 =W (x) dx(4.5.28)aand where vj,1 is the first component of v.

As we shall see in Chapter 11, finding alleigenvalues and eigenvectors of a symmetric tridiagonal matrix is a relatively efficient andwell-conditioned procedure. We accordingly give a routine, gaucof, for finding the abscissasand weights, given the coefficients aj and bj . Remember that if you know the recurrencerelation for orthogonal polynomials that are not normalized to be monic, you can easilyconvert it to monic form by means of the quantities λj .#include <math.h>#include "nrutil.h"void gaucof(int n, float a[], float b[], float amu0, float x[], float w[])Computes the abscissas and weights for a Gaussian quadrature formula from the Jacobi matrix.On input, a[1..n] and b[1..n] are the coefficients of the recurrence relation for the set of'monic orthogonal polynomials.

The quantity µ0 ≡ ab W (x) dx is input as amu0. The abscissasx[1..n] are returned in descending order, with the corresponding weights in w[1..n]. Thearrays a and b are modified. Execution can be speeded up by modifying tqli and eigsrt tocompute only the first component of each eigenvector.{void eigsrt(float d[], float **v, int n);void tqli(float d[], float e[], int n, float **z);int i,j;float **z;z=matrix(1,n,1,n);for (i=1;i<=n;i++) {if (i != 1) b[i]=sqrt(b[i]);Set up superdiagonal of Jacobi matrix.for (j=1;j<=n;j++) z[i][j]=(float)(i == j);Set up identity matrix for tqli to compute eigenvectors.}tqli(a,b,n,z);eigsrt(a,z,n);Sort eigenvalues into descending order.for (i=1;i<=n;i++) {x[i]=a[i];w[i]=amu0*z[1][i]*z[1][i];Equation (4.5.12).}free_matrix(z,1,n,1,n);}Orthogonal Polynomials with Nonclassical WeightsThis somewhat specialized subsection will tell you what to do if your weight functionis not one of the classical ones dealt with above and you do not know the aj ’s and bj ’sof the recurrence relation (4.5.6) to use in gaucof.

Then, a method of finding the aj ’sand bj ’s is needed.The procedure of Stieltjes is to compute a0 from (4.5.7), then p1 (x) from (4.5.6).Knowing p0 and p1 , we can compute a1 and b1 from (4.5.7), and so on. But how are weto compute the inner products in (4.5.7)?158Chapter 4.Integration of FunctionsThe textbook approach is to represent each pj (x) explicitly as a polynomial in x andto compute the inner products by multiplying out term by term. This will be feasible if weknow the first 2N moments of the weight function,& bµj =xj W (x)dxj = 0, 1, . .

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

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