/******************************************************* * Collection of Chebyshev approximation routines * * ---------------------------------------------------- * * REFERENCE: "Numerical Recipes, The Art of Scientific * * Computing By W.H. Press, B.P. Flannery, * * S.A. Teukolsky and W.T. Vetterling, * * Cambridge University Press, 1986" * * [BIBLI 08]. * * * * C++ Release 1.1 By J-P Moreau, Paris. * * (www.jpmoreau.fr) * * ---------------------------------------------------- * * Release 1.1: added functions CHINT and CHDER. * *******************************************************/ #include #include #define NMAX 51 #define HALF 0.5 #define ONE 1.0 #define QUART 0.25 #define TWO 2.0 #define ZERO 0.0 #define PI 4.0*atan(1.0) double FUNC(double x); void CHEBFT(double A, double B, double *C, int N) { /******************************************************* * Chebyshev fit: Given a real function FUNC(X), lower * * and upper limits of the interval [A,B] for X, and a * * maximum degree N, this routine computes the N Cheby- * * shev coefficients Ck, such that FUNC(X) is approxima-* * ted by: N * * [Sum Ck Tk-1(Y)] - C1/2, where X and Y are * * k=1 * * related by: Y = (X - 1/2(A+B)) / (1/2(B-A)) * * This routine is to be used with moderately large N * * (e.g. 30 or 50), the array of C's subsequently to be * * truncated at the smaller value m such that Cm+1 and * * subsequent elements are negligible. * *******************************************************/ double SUM, F[NMAX]; double BMA,BPA,FAC, Y; int J,K; BMA=HALF*(B-A); BPA=HALF*(B+A); for (K=1; K<=N; K++) { Y=cos(PI*(K-HALF)/N); F[K]=FUNC(Y*BMA+BPA); } FAC=TWO/N; for (J=1; J<=N; J++) { SUM=ZERO; for (K=1; K<=N; K++) SUM += F[K]*cos((PI*(J-1))*((K-HALF)/N)); C[J]=FAC*SUM; } } double CHEBEV(double A, double B, double *C, int M, double X) { /********************************************************* * Chebyshev evaluation: All arguments are input. C is an * * array of Chebyshev coefficients, of length M, the first* * M elements of Coutput from subroutine CHEBFT (which * * must have been called with the same A and B). The Che- * * byshev polynomial is evaluated at a point Y determined * * from X, A and B, and the result FUNC(X) is returned as * * the function value. * *********************************************************/ double D,DD,SV,Y,Y2; int J; if ((X-A)*(X-B) > ZERO) printf("\n X not in range.\n\n"); D=ZERO; DD=ZERO; Y=(TWO*X-A-B)/(B-A); // change of variable Y2=TWO*Y; for (J=M; J>1; J--) { SV=D; D=Y2*D-DD+C[J]; DD=SV; } return (Y*D-DD+HALF*C[1]); } void CHINT(double A,double B, double *C, double *CINT, int N) { /********************************************************* * Given A,B,C, as output from routine CHEBFT, and given * * N, the desired degree of approximation (length of C to * * be used), this routine returns the array CINT, the Che-* * byshev coefficients of the integral of the function * * whose coefficients are C. The constant of integration * * is set so that the integral vanishes at A. * *********************************************************/ double CON,FAC,SUM; int J; CON=QUART*(B-A); SUM=ZERO; FAC=ONE; for (J=2; J= 3) for (J=N-2; J>0; J--) CDER[J]=CDER[J+2]+TWO*J*C[J+1]; CON=TWO/(B-A); for (J=1; J<=N; J++) //normalize to interval B - A CDER[J] *= CON; } // end of file Chebyshe.cpp