!******************************************************** !* 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]. * !* * !* F90 Release 1.1 By J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !* ---------------------------------------------------- * !* Release 1.1: added subroutines CHINT and CHDER. * !******************************************************** subroutine CHEBFT(A,B,C,N) real*8 A,B,C(N) external FUNC !******************************************************** !* 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. * !******************************************************** parameter(NMAX=50,HALF=0.5d0,TWO=2.d0,ZERO=0.d0) real*8 PI,SUM,F(NMAX) real*8 BMA,BPA,FAC, Y real*8 FUNC PI=4.d0*datan(1.d0) BMA=HALF*(B-A); BPA=HALF*(B+A) do K=1,N Y=DCOS(PI*(K-HALF)/N) F(K)=FUNC(Y*BMA+BPA) end do FAC=TWO/N do J=1,N SUM=ZERO do K=1,N SUM=SUM+F(K)*DCOS((PI*(J-1))*((K-HALF)/N)) end do C(J)=FAC*SUM end do return end real*8 function CHEBEV(A,B,C,M,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. * !********************************************************** parameter(HALF=0.5d0,TWO=2.d0,ZERO=0.d0) real*8 A,B,C(M),X real*8 D,DD,SV,Y,Y2 if ((X-A)*(X-B).gt.ZERO) pause ' X not in range.' D=ZERO; DD=ZERO Y=(TWO*X-A-B)/(B-A) !change of variable Y2=TWO*Y do J=M,2,-1 SV=D D=Y2*D-DD+C(J) DD=SV end do CHEBEV=Y*D-DD+HALF*C(1) return end subroutine CHINT(A,B,C,CINT,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. * !********************************************************** real*8 A,B,C(N),CINT(N) parameter(ONE=1.d0,QUART=0.25d0,TWO=2.d0,ZERO=0.d0) real*8 CON,FAC,SUM CON=QUART*(B-A) SUM=ZERO FAC=ONE do J=2, N-1 CINT(J)=CON*(C(J-1)-C(J+1))/(J-1) SUM=SUM+FAC*CINT(J) FAC=-FAC end do CINT(N)=CON*C(N-1)/(N-1) SUM=SUM+FAC*CINT(N) CINT(1)=TWO*SUM !set the constant of integration return end subroutine CHDER(A,B,C,CDER,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 CDER, the Che-* !* byshev coefficients of the derivative of the function * !* whose coefficients are C. * !********************************************************** real*8 A,B,C(N),CDER(N) parameter(TWO=2.d0,ZERO=0.d0) real*8 CON CDER(N)=ZERO CDER(N-1)=TWO*(N-1)*C(N) if (N.ge.3) then do J=N-2, 1, -1 CDER(J)=CDER(J+2)+TWO*J*C(J+1) end do end if CON=TWO/(B-A) do J=1,N !normalize to interval B - A CDER(J)=CDER(J)*CON end do return end !end of file Chebyshe.f90