Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 50
Текст из файла (страница 50)
This is effected by a change of variabley≡x − 12 (b + a)12 (b − a)(5.8.10)and by the approximation of f(x) by a Chebyshev polynomial in y.#include <math.h>#include "nrutil.h"#define PI 3.141592653589793void chebft(float a, float b, float c[], int n, float (*func)(float))Chebyshev fit: Given a function func, lower and upper limits of the interval [a,b], and amaximum degree n, this routine computes the n coefficients c[0..n-1] such that func(x) ≈"[ n-1k=0 ck Tk (y)] − c0 /2, where y and x are related by (5.8.10). This routine is to be used withmoderately large n (e.g., 30 or 50), the array of c’s subsequently to be truncated at the smallervalue m such that cm and subsequent elements are negligible.{int k,j;float fac,bpa,bma,*f;5.8 Chebyshev Approximation193f=vector(0,n-1);bma=0.5*(b-a);bpa=0.5*(b+a);for (k=0;k<n;k++) {We evaluate the function at the n points requiredfloat y=cos(PI*(k+0.5)/n);by (5.8.7).f[k]=(*func)(y*bma+bpa);}fac=2.0/n;for (j=0;j<n;j++) {double sum=0.0;We will accumulate the sum in double precision,for (k=0;k<n;k++)a nicety that you can ignore.sum += f[k]*cos(PI*j*(k+0.5)/n);c[j]=fac*sum;}free_vector(f,0,n-1);}(If you find that the execution time of chebft is dominated by the calculation ofN 2 cosines, rather than by the N evaluations of your function, then you should lookahead to §12.3, especially equation 12.3.22, which shows how fast cosine transformmethods can be used to evaluate equation 5.8.7.)Now that we have the Chebyshev coefficients, how do we evaluate the approximation? One could use the recurrence relation of equation (5.8.2) to generate valuesfor Tk (x) from T0 = 1, T1 = x, while also accumulating the sum of (5.8.9).
Itis better to use Clenshaw’s recurrence formula (§5.5), effecting the two processessimultaneously. Applied to the Chebyshev series (5.8.9), the recurrence isdm+1 ≡ dm ≡ 0dj = 2xdj+1 − dj+2 + cjj = m − 1, m − 2, . . . , 1(5.8.11)1f(x) ≡ d0 = xd1 − d2 + c02float chebev(float a, float b, float c[], int m, float x)Chebyshev evaluation: All arguments are input. c[0..m-1] is an array of Chebyshev coefficients, the first m elements of c output from chebft (which must have been called with the"m-1same a and b). The Chebyshev polynomialk=0 ck Tk (y) − c0 /2 is evaluated at a pointy = [x − (b + a)/2]/[(b − a)/2], and the result is returned as the function value.{void nrerror(char error_text[]);float d=0.0,dd=0.0,sv,y,y2;int j;if ((x-a)*(x-b) > 0.0) nrerror("x not in range in routine chebev");y2=2.0*(y=(2.0*x-a-b)/(b-a));Change of variable.for (j=m-1;j>=1;j--) {Clenshaw’s recurrence.sv=d;d=y2*d-dd+c[j];dd=sv;}return y*d-dd+0.5*c[0];Last step is different.}194Chapter 5.Evaluation of FunctionsIf we are approximating an even function on the interval [−1, 1], its expansionwill involve only even Chebyshev polynomials.
It is wasteful to call chebev withall the odd coefficients zero [1]. Instead, using the half-angle identity for the cosinein equation (5.8.1), we get the relationT2n (x) = Tn (2x2 − 1)(5.8.12)Thus we can evaluate a series of even Chebyshev polynomials by calling chebevwith the even coefficients stored consecutively in the array c, but with the argumentx replaced by 2x2 − 1.An odd function will have an expansion involving only odd Chebysev polynomials.
It is best to rewrite it as an expansion for the function f(x)/x, whichinvolves only even Chebyshev polynomials. This will give accurate values forf(x)/x near x = 0. The coefficients cn for f(x)/x can be found from those forf(x) by recurrence:cN+1 = 0cn−1 = 2cn − cn+1 ,n = N, N − 2, . . .(5.8.13)Equation (5.8.13) follows from the recurrence relation in equation (5.8.2).If you insist on evaluating an odd Chebyshev series, the efficient way is to onceagain use chebev with x replaced by y = 2x2 − 1, and with the odd coefficientsstored consecutively in the array c. Now, however, you must also change the lastformula in equation (5.8.11) to bef(x) = x[(2y − 1)d2 − d3 + c1 ](5.8.14)and change the corresponding line in chebev.CITED REFERENCES AND FURTHER READING:Clenshaw, C.W. 1962, Mathematical Tables, vol. 5, National Physical Laboratory, (London: H.M.Stationery Office).
[1]Goodwin, E.T. (ed.) 1961, Modern Computing Methods, 2nd ed. (New York: Philosophical Library), Chapter 8.Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall),§4.4.1, p. 104.Johnson, L.W., and Riess, R.D. 1982, Numerical Analysis, 2nd ed. (Reading, MA: AddisonWesley), §6.5.2, p. 334.Carnahan, B., Luther, H.A., and Wilkes, J.O. 1969, Applied Numerical Methods (New York:Wiley), §1.10, p. 39.5.9 Derivatives or Integrals of a Chebyshev-approximated Function1955.9 Derivatives or Integrals of aChebyshev-approximated FunctionIf you have obtained the Chebyshev coefficients that approximate a function ina certain range (e.g., from chebft in §5.8), then it is a simple matter to transformthem to Chebyshev coefficients corresponding to the derivative or integral of thefunction.
Having done this, you can evaluate the derivative or integral just as if itwere a function that you had Chebyshev-fitted ab initio.The relevant formulas are these: If ci , i = 0, . . . , m − 1 are the coefficientsthat approximate a function f in equation (5.8.9), Ci are the coefficients thatapproximate the indefinite integral of f, and ci are the coefficients that approximatethe derivative of f, thenCi =ci−1 − ci+12(i − 1)ci−1 = ci+1 + 2(i − 1)ci(i > 1)(5.9.1)(i = m − 1, m − 2, .
. . , 2)(5.9.2)Equation (5.9.1) is augmented by an arbitrary choice of C0 , corresponding to anarbitrary constant of integration. Equation (5.9.2), which is a recurrence, is startedwith the values cm = cm−1 = 0, corresponding to no information about the m + 1stChebyshev coefficient of the original function f.Here are routines for implementing equations (5.9.1) and (5.9.2).void chder(float a, float b, float c[], float cder[], int n)Given a,b,c[0..n-1], as output from routine chebft §5.8, and given n, the desired degreeof approximation (length of c to be used), this routine returns the array cder[0..n-1], theChebyshev coefficients of the derivative of the function whose coefficients are c.{int j;float con;cder[n-1]=0.0;cder[n-2]=2*(n-1)*c[n-1];for (j=n-3;j>=0;j--)cder[j]=cder[j+2]+2*(j+1)*c[j+1];con=2.0/(b-a);for (j=0;j<n;j++)cder[j] *= con;n-1 and n-2 are special cases.Equation (5.9.2).Normalize to the interval b-a.}void chint(float a, float b, float c[], float cint[], int n)Given a,b,c[0..n-1], as output from routine chebft §5.8, and given n, the desired degreeof approximation (length of c to be used), this routine returns the array cint[0..n-1], theChebyshev coefficients of the integral of the function whose coefficients are c.
The constant ofintegration is set so that the integral vanishes at a.{int j;float sum=0.0,fac=1.0,con;con=0.25*(b-a);for (j=1;j<=n-2;j++) {cint[j]=con*(c[j-1]-c[j+1])/j;Factor that normalizes to the interval b-a.Equation (5.9.1).196Chapter 5.Evaluation of Functionssum += fac*cint[j];fac = -fac;}cint[n-1]=con*c[n-2]/(n-1);sum += fac*cint[n-1];cint[0]=2.0*sum;Accumulates the constant of integration.Will equal ±1.Special case of (5.9.1) for n-1.Set the constant of integration.}Clenshaw-Curtis QuadratureSince a smooth function’s Chebyshev coefficients ci decrease rapidly, generally exponentially, equation (5.9.1) is often quite efficient as the basis for a quadrature scheme. Theroutines' x chebft and chint, used in that order, can be followed by repeated calls to chebevif a f (x)dx is required for many different values of x in the range a ≤ x ≤ b.'bIf only the single definite integral a f (x)dx is required, then chint and chebev arereplaced by the simpler formula, derived from equation (5.9.1),& b1111f (x)dx = (b − a)c1 − c3 −c5 − · · · −c2k+1 − · · ·2315(2k + 1)(2k − 1)a(5.9.3)where the ci ’s are as returned by chebft.
The series can be truncated when c2k+1 becomesnegligible, and the first neglected term gives an error estimate.This scheme is known as Clenshaw-Curtis quadrature [1]. It is often combined with anadaptive choice of N , the number of Chebyshev coefficients calculated via equation (5.8.7),which is also the number of function evaluations of f (x). If a modest choice of N doesnot give a sufficiently small c2k+1 in equation (5.9.3), then a larger value is tried.
In thisadaptive case, it is even better to replace equation (5.8.7) by the so-called “trapezoidal” orGauss-Lobatto (§4.5) variant,cj = Nπkπ(j − 1)k2 f coscosNNNj = 1, . . . , N(5.9.4)k=0where (N.B.!) the two primes signify that the first and last terms in the sum are to bemultiplied by 1/2. If N is doubled in equation (5.9.4), then half of the new functionevaluation points are identical to the old ones, allowing the previous function evaluations to bereused. This feature, plus the analytic weights and abscissas (cosine functions in 5.9.4), giveClenshaw-Curtis quadrature an edge over high-order adaptive Gaussian quadrature (cf. §4.5),which the method otherwise resembles.If your problem forces you to large values of N , you should be aware that equation (5.9.4)can be evaluated rapidly, and simultaneously for all the values of j, by a fast cosine transform.(See §12.3, especially equation 12.3.17.) (We already remarked that the nontrapezoidal form(5.8.7) can also be done by fast cosine methods, cf.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
