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

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

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

equation 12.3.22.)CITED REFERENCES AND FURTHER READING:Goodwin, E.T. (ed.) 1961, Modern Computing Methods, 2nd ed. (New York: Philosophical Library), pp. 78–79.Clenshaw, C.W., and Curtis, A.R. 1960, Numerische Mathematik, vol. 2, pp. 197–205. [1]5.10 Polynomial Approximation from Chebyshev Coefficients1975.10 Polynomial Approximation fromChebyshev CoefficientsYou may well ask after reading the preceding two sections, “Must I store andevaluate my Chebyshev approximation as an array of Chebyshev coefficients for atransformed variable y? Can’t I convert the ck ’s into actual polynomial coefficientsin the original variable x and have an approximation of the following form?”f(x) ≈m−1gk xk(5.10.1)k=0Yes, you can do this (and we will give you the algorithm to do it), but wecaution you against it: Evaluating equation (5.10.1), where the coefficient g’s reflectan underlying Chebyshev approximation, usually requires more significant figuresthan evaluation of the Chebyshev sum directly (as by chebev).

This is becausethe Chebyshev polynomials themselves exhibit a rather delicate cancellation: Theleading coefficient of Tn (x), for example, is 2n−1 ; other coefficients of Tn (x) areeven bigger; yet they all manage to combine into a polynomial that lies between ±1.Only when m is no larger than 7 or 8 should you contemplate writing a Chebyshevfit as a direct polynomial, and even in those cases you should be willing to toleratetwo or so significant figures less accuracy than the roundoff limit of your machine.You get the g’s in equation (5.10.1) from the c’s output from chebft (suitablytruncated at a modest value of m) by calling in sequence the following two procedures:#include "nrutil.h"void chebpc(float c[], float d[], int n)Chebyshev polynomial coefficients.

Given a coefficient array c[0..n-1], this routine generates""n-1ka coefficient array d[0..n-1] such that n-1k=0 dk y =k=0 ck Tk (y) − c0 /2. The methodis Clenshaw’s recurrence (5.8.11), but now applied algebraically rather than arithmetically.{int k,j;float sv,*dd;dd=vector(0,n-1);for (j=0;j<n;j++) d[j]=dd[j]=0.0;d[0]=c[n-1];for (j=n-2;j>=1;j--) {for (k=n-j;k>=1;k--) {sv=d[k];d[k]=2.0*d[k-1]-dd[k];dd[k]=sv;}sv=d[0];d[0] = -dd[0]+c[j];dd[0]=sv;}for (j=n-1;j>=1;j--)d[j]=d[j-1]-dd[j];d[0] = -dd[0]+0.5*c[0];free_vector(dd,0,n-1);}198Chapter 5.Evaluation of Functionsvoid pcshft(float a, float b, float d[], int n)Polynomial coefficient shift. Given a coefficient array d[0..n-1], this routine generates a""n-1kkcoefficient array g [0..n-1] such that n-1k=0 dk y =k=0 gk x , where x and y are relatedby (5.8.10), i.e., the interval −1 < y < 1 is mapped to the interval a < x < b.

The arrayg is returned in d.{int k,j;float fac,cnst;cnst=2.0/(b-a);fac=cnst;for (j=1;j<n;j++) {First we rescale by the factor const...d[j] *= fac;fac *= cnst;}cnst=0.5*(a+b);...which is then redefined as the desired shift.for (j=0;j<=n-2;j++)We accomplish the shift by synthetic division. Syntheticfor (k=n-2;k>=j;k--)division is a miracle of high-school algebra.

If youd[k] -= cnst*d[k+1];never learned it, go do so. You won’t be sorry.}CITED REFERENCES AND FURTHER READING:Acton, F.S. 1970, Numerical Methods That Work; 1990, corrected edition (Washington: Mathematical Association of America), pp. 59, 182–183 [synthetic division].5.11 Economization of Power SeriesOne particular application of Chebyshev methods, the economization of power series, isan occasionally useful technique, with a flavor of getting something for nothing.Suppose that you are already computing a function by the use of a convergent powerseries, for examplef (x) ≡ 1 −xx2x3+−+···3!5!7!(5.11.1)√ √(This function is actually sin( x)/ x, but pretend you don’t know that.) You might bedoing a problem that requires evaluating the series many times in some particular interval, say[0, (2π)2 ].

Everything is fine, except that the series requires a large number of terms beforeits error (approximated by the first neglected term, say) is tolerable. In our example, withx = (2π)2 , the first term smaller than 10−7 is x13 /(27!). This then approximates the errorof the finite series whose last term is x12 /(25!).Notice that because of the large exponent in x13 , the error is much smaller than 10−7everywhere in the interval except at the very largest values of x.

This is the feature that allows“economization”: if we are willing to let the error elsewhere in the interval rise to about thesame value that the first neglected term has at the extreme end of the interval, then we canreplace the 13-term series by one that is significantly shorter.Here are the steps for doing so:1. Change variables from x to y, as in equation (5.8.10), to map the x interval into−1 ≤ y ≤ 1.2. Find the coefficients of the Chebyshev sum (like equation 5.8.8) that exactly equals yourtruncated power series (the one with enough terms for accuracy).3. Truncate this Chebyshev series to a smaller number of terms, using the coefficient of thefirst neglected Chebyshev polynomial as an estimate of the error.5.11 Economization of Power Series1994. Convert back to a polynomial in y.5.

Change variables back to x.All of these steps can be done numerically, given the coefficients of the original powerseries expansion. The first step is exactly the inverse of the routine pcshft (§5.10), whichmapped a polynomial from y (in the interval [−1, 1]) to x (in the interval [a, b]). But sinceequation (5.8.10) is a linear relation between x and y, one can also use pcshft for theinverse. The inverse ofpcshft(a,b,d,n)turns out to be (you can check this)−2 − b − a 2 − b − apcshft,,d,nb−ab−aThe second step requires the inverse operation to that done by the routine chebpc (whichtook Chebyshev coefficients into polynomial coefficients).

The following routine, pccheb,accomplishes this, using the formula [1] !1kkxk = k−1 Tk (x) +Tk−2 (x) +Tk−4 (x) + · · ·(5.11.2)122where the last term depends on whether k is even or odd,k1 k··· +T1 (x) (k odd),T0 (x) (k even).··· +(k − 1)/22 k/2(5.11.3)void pccheb(float d[], float c[], int n)Inverse of routine chebpc: given an array of polynomial coefficients d[0..n-1], returns anequivalent array of Chebyshev coefficients c[0..n-1].{int j,jm,jp,k;float fac,pow;pow=1.0;Will be powers of 2.c[0]=2.0*d[0];for (k=1;k<n;k++) {Loop over orders of x in the polynomial.c[k]=0.0;Zero corresponding order of Chebyshev.fac=d[k]/pow;jm=k;jp=1;for (j=k;j>=0;j-=2,jm--,jp++) {Increment this and lower orders of Chebyshev with the combinatorial coefficent timesd[k]; see text for formula.c[j] += fac;fac *= ((float)jm)/((float)jp);}pow += pow;}}The fourth and fifth steps are accomplished by the routines chebpc and pcshft,respectively.

Here is how the procedure looks all together:200Chapter 5.Evaluation of Functions#define NFEW ..#define NMANY ..float *c,*d,*e,a,b;Economize NMANY power series coefficients e[0..NMANY-1] in the range (a, b) into NFEWcoefficients d[0..NFEW-1].c=vector(0,NMANY-1);d=vector(0,NFEW-1);e=vector(0,NMANY-1);pcshft((-2.0-b-a)/(b-a),(2.0-b-a)/(b-a),e,NMANY);pccheb(e,c,NMANY);...Here one would normally examine the Chebyshev coefficients c[0..NMANY-1] to decidehow small NFEW can be.chebpc(c,d,NFEW);pcshft(a,b,d,NFEW);In our example, by the way, the 8th through 10th Chebyshev coefficients turn out tobe on the order of −7 × 10−6 , 3 × 10−7 , and −9 × 10−9 , so reasonable truncations (forsingle precision calculations) are somewhere in this range, yielding a polynomial with 8 –10 terms instead of the original 13.Replacing a 13-term polynomial with a (say) 10-term polynomial without any loss ofaccuracy — that does seem to be getting something for nothing.

Is there some magic inthis technique? Not really. The 13-term polynomial defined a function f (x). Equivalent toeconomizing the series, we could instead have evaluated f (x) at enough points to constructits Chebyshev approximation in the interval of interest, by the methods of §5.8. We wouldhave obtained just the same lower-order polynomial. The principal lesson is that the rateof convergence of Chebyshev coefficients has nothing to do with the rate of convergence ofpower series coefficients; and it is the former that dictates the number of terms needed in apolynomial approximation. A function might have a divergent power series in some regionof interest, but if the function itself is well-behaved, it will have perfectly good polynomialapproximations.

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

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