c4-5 (779482), страница 4
Текст из файла (страница 4)
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)?Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).where(4.5.27)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,Z bµj =xj W (x)dxj = 0, 1, . . . , 2N − 1(4.5.29)aawhere the πj ’s satisfy a recurrence relation analogous to (4.5.6),π−1 (x) ≡ 0π0 (x) ≡ 1(4.5.31)πj+1 (x) = (x − αj )πj (x) − βj πj−1 (x)j = 0, 1, 2, .
. .and the coefficients αj , βj are known explicitly. Then Wheeler [6] has given an efficientO(N 2 ) algorithm equivalent to that of Sack and Donovan for finding aj and bj via a setof intermediate quantitiesσk,l = hpk |πl ik, l ≥ −1(4.5.32)Initializel = 1, 2, . . . , 2N − 2l = 0, 1, .
. . , 2N − 1σ−1,l = 0σ0,l = νla0 = α0 +(4.5.33)ν1ν0b0 = 0Then, for k = 1, 2, . . . , N − 1, computeσk,l = σk−1,l+1 − (ak−1 − αl )σk−1,l − bk−1 σk−2,l + βlσk−1,l−1l = k, k + 1, . . . , 2N − k − 1σk−1,kσk,k+1ak = αk −+σk−1,k−1σk,kσk,kbk =σk−1,k−1(4.5.34)Note that the normalization factors can also easily be computed if needed:hp0 |p0 i = ν0hpj |pj i = bj hpj−1 |pj−1 ij = 1, 2, . . .(4.5.35)You can find a derivation of the above algorithm in Ref. [7].Wheeler’s algorithm requires that the modified moments (4.5.30) be accurately computed.In practical cases there is often a closed form, or else recurrence relations can be used. TheSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).However, the solution of the resulting set of algebraic equations for the coefficients aj and bjin terms of the moments µj is in general extremely ill-conditioned. Even in double precision,it is not unusual to lose all accuracy by the time N = 12. We thus reject any procedurebased on the moments (4.5.29).Sack and Donovan [5] discovered that the numerical stability is greatly improved if,instead of using powers of x as a set of basis functions to represent the pj ’s, one uses someother known set of orthogonal polynomials πj (x), say. Roughly speaking, the improvedstability occurs because the polynomial basis “samples” the interval (a, b) better than thepower basis when the inner product integrals are evaluated, especially if its weight functionresembles W (x).So assume that we know the modified momentsZ bνj =πj (x)W (x)dxj = 0, 1, .
. . , 2N − 1(4.5.30)4.5 Gaussian Quadratures and Orthogonal Polynomials159#include "nrutil.h"void orthog(int n, float anu[], float alpha[], float beta[], float a[],float b[])Computes the coefficients aj and bj , j = 0, . . . N − 1, of the recurrence relation for monicorthogonal polynomials with weight function W (x) by Wheeler’s algorithm. On input, the arraysalpha[1..2*n-1] and beta[1..2*n-1] are the coefficients αj and βj , j = 0, .
. . 2N −2, ofthe recurrence relation for the chosen basis of orthogonal polynomials. The modified momentsνj are input in anu[1..2*n]. The first n coefficients are returned in a[1..n] and b[1..n].{int k,l;float **sig;int looptmp;sig=matrix(1,2*n+1,1,2*n+1);looptmp=2*n;for (l=3;l<=looptmp;l++) sig[1][l]=0.0;Initialization, Equation (4.5.33).looptmp++;for (l=2;l<=looptmp;l++) sig[2][l]=anu[l-1];a[1]=alpha[1]+anu[2]/anu[1];b[1]=0.0;for (k=3;k<=n+1;k++) {Equation (4.5.34).looptmp=2*n-k+3;for (l=k;l<=looptmp;l++) {sig[k][l]=sig[k-1][l+1]+(alpha[l-1]-a[k-2])*sig[k-1][l]b[k-2]*sig[k-2][l]+beta[l-1]*sig[k-1][l-1];}a[k-1]=alpha[k-1]+sig[k][k+1]/sig[k][k]-sig[k-1][k]/sig[k-1][k-1];b[k-1]=sig[k][k]/sig[k-1][k-1];}free_matrix(sig,1,2*n+1,1,2*n+1);}As an example of the use of orthog, consider the problem [7] of generating orthogonalpolynomials with the weight function W (x) = − log x on the interval (0, 1).
A suitable setof πj ’s is the shifted Legendre polynomialsπj =(j!)2Pj (2x − 1)(2j)!(4.5.36)The factor in front of Pj makes the polynomials monic. The coefficients in the recurrencerelation (4.5.31) are1αj =j = 0, 1, . . .2(4.5.37)1βj =j=1,2,...4(4 − j −2 )while the modified moments areνj =1(−1)j (j!)2j(j + 1)(2j)!j =0j ≥1(4.5.38)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).algorithm is extremely successful for finite intervals (a, b). For infinite intervals, the algorithmdoes not completely remove the ill-conditioning. In this case, Gautschi [8,9] recommendsreducing the interval to a finite interval by a change of variable, and then using a suitablediscretization procedure to compute the inner products. You will have to consult thereferences for details.We give the routine orthog for generating the coefficients aj and bj by Wheeler’salgorithm, given the coefficients αj and βj , and the modified moments νj . For consistencywith gaucof, the vectors α, β, a and b are 1-based.
Correspondingly, we increase the indicesof the σ matrix by 2, i.e., sig[k,l] = σk−2,l−2 .160Chapter 4.Integration of FunctionsA call to orthog with this input allows one to generate the required polynomials to machineaccuracy for very large N , and hence do Gaussian quadrature with this weight function.
BeforeSack and Donovan’s observation, this seemingly simple problem was essentially intractable.Extensions of Gaussian QuadratureCITED REFERENCES AND FURTHER READING:Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions, Applied Mathematics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover Publications, New York), §25.4. [1]Stroud, A.H., and Secrest, D. 1966, Gaussian Quadrature Formulas (Englewood Cliffs, NJ:Prentice-Hall). [2]Golub, G.H., and Welsch, J.H.
1969, Mathematics of Computation, vol. 23, pp. 221–230 andA1–A10. [3]Wilf, H.S. 1962, Mathematics for the Physical Sciences (New York: Wiley), Problem 9, p. 80. [4]Sack, R.A., and Donovan, A.F. 1971/72, Numerische Mathematik, vol. 18, pp. 465–478. [5]Wheeler, J.C. 1974, Rocky Mountain Journal of Mathematics, vol. 4, pp. 287–296. [6]Gautschi, W. 1978, in Recent Advances in Numerical Analysis, C. de Boor and G.H. Golub, eds.(New York: Academic Press), pp.
45–72. [7]Gautschi, W. 1981, in E.B. Christoffel, P.L. Butzer and F. Fehér, eds. (Basel: Birkhauser Verlag),pp. 72–147. [8]Gautschi, W. 1990, in Orthogonal Polynomials, P. Nevai, ed. (Dordrecht: Kluwer AcademicPublishers), pp. 181–216.
[9]Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).There are many different ways in which the ideas of Gaussian quadrature havebeen extended. One important extension is the case of preassigned nodes: Somepoints are required to be included in the set of abscissas, and the problem is to choosethe weights and the remaining abscissas to maximize the degree of exactness of thethe quadrature rule.
The most common cases are Gauss-Radau quadrature, whereone of the nodes is an endpoint of the interval, either a or b, and Gauss-Lobattoquadrature, where both a and b are nodes. Golub [10] has given an algorithm similarto gaucof for these cases.The second important extension is the Gauss-Kronrod formulas. For ordinaryGaussian quadrature formulas, as N increases the sets of abscissas have no pointsin common. This means that if you compare results with increasing N as a way ofestimating the quadrature error, you cannot reuse the previous function evaluations.Kronrod [11] posed the problem of searching for optimal sequences of rules, eachof which reuses all abscissas of its predecessor. If one starts with N = m, say,and then adds n new points, one has 2n + m free parameters: the n new abscissasand weights, and m new weights for the fixed previous abscissas.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
