*DECK DQC25F SUBROUTINE DQC25F (F, A, B, OMEGA, INTEGR, NRMOM, MAXP1, KSAVE, + RESULT, ABSERR, NEVAL, RESABS, RESASC, MOMCOM, CHEBMO) C***BEGIN PROLOGUE DQC25F C***PURPOSE To compute the integral I=Integral of F(X) over (A,B) C Where W(X) = COS(OMEGA*X) or W(X)=SIN(OMEGA*X) and to C compute J = Integral of ABS(F) over (A,B). For small value C of OMEGA or small intervals (A,B) the 15-point GAUSS-KRONRO C Rule is used. Otherwise a generalized CLENSHAW-CURTIS C method is used. C***LIBRARY SLATEC (QUADPACK) C***CATEGORY H2A2A2 C***TYPE DOUBLE PRECISION (QC25F-S, DQC25F-D) C***KEYWORDS CLENSHAW-CURTIS METHOD, GAUSS-KRONROD RULES, C INTEGRATION RULES FOR FUNCTIONS WITH COS OR SIN FACTOR, C QUADPACK, QUADRATURE C***AUTHOR Piessens, Robert C Applied Mathematics and Programming Division C K. U. Leuven C de Doncker, Elise C Applied Mathematics and Programming Division C K. U. Leuven C***DESCRIPTION C C Integration rules for functions with COS or SIN factor C Standard fortran subroutine C Double precision version C C PARAMETERS C ON ENTRY C F - Double precision C Function subprogram defining the integrand C function F(X). The actual name for F needs to C be declared E X T E R N A L in the calling program. C C A - Double precision C Lower limit of integration C C B - Double precision C Upper limit of integration C C OMEGA - Double precision C Parameter in the WEIGHT function C C INTEGR - Integer C Indicates which WEIGHT function is to be used C INTEGR = 1 W(X) = COS(OMEGA*X) C INTEGR = 2 W(X) = SIN(OMEGA*X) C C NRMOM - Integer C The length of interval (A,B) is equal to the length C of the original integration interval divided by C 2**NRMOM (we suppose that the routine is used in an C adaptive integration process, otherwise set C NRMOM = 0). NRMOM must be zero at the first call. C C MAXP1 - Integer C Gives an upper bound on the number of Chebyshev C moments which can be stored, i.e. for the C intervals of lengths ABS(BB-AA)*2**(-L), C L = 0,1,2, ..., MAXP1-2. C C KSAVE - Integer C Key which is one when the moments for the C current interval have been computed C C ON RETURN C RESULT - Double precision C Approximation to the integral I C C ABSERR - Double precision C Estimate of the modulus of the absolute C error, which should equal or exceed ABS(I-RESULT) C C NEVAL - Integer C Number of integrand evaluations C C RESABS - Double precision C Approximation to the integral J C C RESASC - Double precision C Approximation to the integral of ABS(F-I/(B-A)) C C ON ENTRY AND RETURN C MOMCOM - Integer C For each interval length we need to compute the C Chebyshev moments. MOMCOM counts the number of C intervals for which these moments have already been C computed. If NRMOM.LT.MOMCOM or KSAVE = 1, the C Chebyshev moments for the interval (A,B) have C already been computed and stored, otherwise we C compute them and we increase MOMCOM. C C CHEBMO - Double precision C Array of dimension at least (MAXP1,25) containing C the modified Chebyshev moments for the first MOMCOM C MOMCOM interval lengths C C ...................................................................... C C***REFERENCES (NONE) C***ROUTINES CALLED D1MACH, DGTSL, DQCHEB, DQK15W, DQWGTF C***REVISION HISTORY (YYMMDD) C 810101 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890531 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C***END PROLOGUE DQC25F C DOUBLE PRECISION A,ABSERR,AC,AN,AN2,AS,ASAP,ASS,B,CENTR,CHEBMO, 1 CHEB12,CHEB24,CONC,CONS,COSPAR,D,DQWGTF,D1, 2 D1MACH,D2,ESTC,ESTS,F,FVAL,HLGTH,OFLOW,OMEGA,PARINT,PAR2,PAR22, 3 P2,P3,P4,RESABS,RESASC,RESC12,RESC24,RESS12,RESS24,RESULT, 4 SINPAR,V,X INTEGER I,IERS,INTEGR,ISYM,J,K,KSAVE,M,MOMCOM,NEVAL,MAXP1, 1 NOEQU,NOEQ1,NRMOM C DIMENSION CHEBMO(MAXP1,25),CHEB12(13),CHEB24(25),D(25),D1(25), 1 D2(25),FVAL(25),V(28),X(11) C EXTERNAL F, DQWGTF C C THE VECTOR X CONTAINS THE VALUES COS(K*PI/24) C K = 1, ...,11, TO BE USED FOR THE CHEBYSHEV EXPANSION OF F C SAVE X DATA X(1),X(2),X(3),X(4),X(5),X(6),X(7),X(8),X(9),X(10),X(11)/ 1 0.9914448613738104D+00, 0.9659258262890683D+00, 2 0.9238795325112868D+00, 0.8660254037844386D+00, 3 0.7933533402912352D+00, 0.7071067811865475D+00, 4 0.6087614290087206D+00, 0.5000000000000000D+00, 5 0.3826834323650898D+00, 0.2588190451025208D+00, 6 0.1305261922200516D+00/ C C LIST OF MAJOR VARIABLES C ----------------------- C C CENTR - MID POINT OF THE INTEGRATION INTERVAL C HLGTH - HALF-LENGTH OF THE INTEGRATION INTERVAL C FVAL - VALUE OF THE FUNCTION F AT THE POINTS C (B-A)*0.5*COS(K*PI/12) + (B+A)*0.5, K = 0, ..., 24 C CHEB12 - COEFFICIENTS OF THE CHEBYSHEV SERIES EXPANSION C OF DEGREE 12, FOR THE FUNCTION F, IN THE C INTERVAL (A,B) C CHEB24 - COEFFICIENTS OF THE CHEBYSHEV SERIES EXPANSION C OF DEGREE 24, FOR THE FUNCTION F, IN THE C INTERVAL (A,B) C RESC12 - APPROXIMATION TO THE INTEGRAL OF C COS(0.5*(B-A)*OMEGA*X)*F(0.5*(B-A)*X+0.5*(B+A)) C OVER (-1,+1), USING THE CHEBYSHEV SERIES C EXPANSION OF DEGREE 12 C RESC24 - APPROXIMATION TO THE SAME INTEGRAL, USING THE C CHEBYSHEV SERIES EXPANSION OF DEGREE 24 C RESS12 - THE ANALOGUE OF RESC12 FOR THE SINE C RESS24 - THE ANALOGUE OF RESC24 FOR THE SINE C C C MACHINE DEPENDENT CONSTANT C -------------------------- C C OFLOW IS THE LARGEST POSITIVE MAGNITUDE. C C***FIRST EXECUTABLE STATEMENT DQC25F OFLOW = D1MACH(2) C CENTR = 0.5D+00*(B+A) HLGTH = 0.5D+00*(B-A) PARINT = OMEGA*HLGTH C C COMPUTE THE INTEGRAL USING THE 15-POINT GAUSS-KRONROD C FORMULA IF THE VALUE OF THE PARAMETER IN THE INTEGRAND C IS SMALL. C IF(ABS(PARINT).GT.0.2D+01) GO TO 10 CALL DQK15W(F,DQWGTF,OMEGA,P2,P3,P4,INTEGR,A,B,RESULT, 1 ABSERR,RESABS,RESASC) NEVAL = 15 GO TO 170 C C COMPUTE THE INTEGRAL USING THE GENERALIZED CLENSHAW- C CURTIS METHOD. C 10 CONC = HLGTH*COS(CENTR*OMEGA) CONS = HLGTH*SIN(CENTR*OMEGA) RESASC = OFLOW NEVAL = 25 C C CHECK WHETHER THE CHEBYSHEV MOMENTS FOR THIS INTERVAL C HAVE ALREADY BEEN COMPUTED. C IF(NRMOM.LT.MOMCOM.OR.KSAVE.EQ.1) GO TO 120 C C COMPUTE A NEW SET OF CHEBYSHEV MOMENTS. C M = MOMCOM+1 PAR2 = PARINT*PARINT PAR22 = PAR2+0.2D+01 SINPAR = SIN(PARINT) COSPAR = COS(PARINT) C C COMPUTE THE CHEBYSHEV MOMENTS WITH RESPECT TO COSINE. C V(1) = 0.2D+01*SINPAR/PARINT V(2) = (0.8D+01*COSPAR+(PAR2+PAR2-0.8D+01)*SINPAR/PARINT)/PAR2 V(3) = (0.32D+02*(PAR2-0.12D+02)*COSPAR+(0.2D+01* 1 ((PAR2-0.80D+02)*PAR2+0.192D+03)*SINPAR)/PARINT)/(PAR2*PAR2) AC = 0.8D+01*COSPAR AS = 0.24D+02*PARINT*SINPAR IF(ABS(PARINT).GT.0.24D+02) GO TO 30 C C COMPUTE THE CHEBYSHEV MOMENTS AS THE SOLUTIONS OF A C BOUNDARY VALUE PROBLEM WITH 1 INITIAL VALUE (V(3)) AND 1 C END VALUE (COMPUTED USING AN ASYMPTOTIC FORMULA). C NOEQU = 25 NOEQ1 = NOEQU-1 AN = 0.6D+01 DO 20 K = 1,NOEQ1 AN2 = AN*AN D(K) = -0.2D+01*(AN2-0.4D+01)*(PAR22-AN2-AN2) D2(K) = (AN-0.1D+01)*(AN-0.2D+01)*PAR2 D1(K+1) = (AN+0.3D+01)*(AN+0.4D+01)*PAR2 V(K+3) = AS-(AN2-0.4D+01)*AC AN = AN+0.2D+01 20 CONTINUE AN2 = AN*AN D(NOEQU) = -0.2D+01*(AN2-0.4D+01)*(PAR22-AN2-AN2) V(NOEQU+3) = AS-(AN2-0.4D+01)*AC V(4) = V(4)-0.56D+02*PAR2*V(3) ASS = PARINT*SINPAR ASAP = (((((0.210D+03*PAR2-0.1D+01)*COSPAR-(0.105D+03*PAR2 1 -0.63D+02)*ASS)/AN2-(0.1D+01-0.15D+02*PAR2)*COSPAR 2 +0.15D+02*ASS)/AN2-COSPAR+0.3D+01*ASS)/AN2-COSPAR)/AN2 V(NOEQU+3) = V(NOEQU+3)-0.2D+01*ASAP*PAR2*(AN-0.1D+01)* 1 (AN-0.2D+01) C C SOLVE THE TRIDIAGONAL SYSTEM BY MEANS OF GAUSSIAN C ELIMINATION WITH PARTIAL PIVOTING. C C *** CALL TO DGTSL MUST BE REPLACED BY CALL TO C *** DOUBLE PRECISION VERSION OF LINPACK ROUTINE SGTSL C CALL DGTSL(NOEQU,D1,D,D2,V(4),IERS) GO TO 50 C C COMPUTE THE CHEBYSHEV MOMENTS BY MEANS OF FORWARD C RECURSION. C 30 AN = 0.4D+01 DO 40 I = 4,13 AN2 = AN*AN V(I) = ((AN2-0.4D+01)*(0.2D+01*(PAR22-AN2-AN2)*V(I-1)-AC) 1 +AS-PAR2*(AN+0.1D+01)*(AN+0.2D+01)*V(I-2))/ 2 (PAR2*(AN-0.1D+01)*(AN-0.2D+01)) AN = AN+0.2D+01 40 CONTINUE 50 DO 60 J = 1,13 CHEBMO(M,2*J-1) = V(J) 60 CONTINUE C C COMPUTE THE CHEBYSHEV MOMENTS WITH RESPECT TO SINE. C V(1) = 0.2D+01*(SINPAR-PARINT*COSPAR)/PAR2 V(2) = (0.18D+02-0.48D+02/PAR2)*SINPAR/PAR2 1 +(-0.2D+01+0.48D+02/PAR2)*COSPAR/PARINT AC = -0.24D+02*PARINT*COSPAR AS = -0.8D+01*SINPAR IF(ABS(PARINT).GT.0.24D+02) GO TO 80 C C COMPUTE THE CHEBYSHEV MOMENTS AS THE SOLUTIONS OF A BOUNDARY C VALUE PROBLEM WITH 1 INITIAL VALUE (V(2)) AND 1 END VALUE C (COMPUTED USING AN ASYMPTOTIC FORMULA). C AN = 0.5D+01 DO 70 K = 1,NOEQ1 AN2 = AN*AN D(K) = -0.2D+01*(AN2-0.4D+01)*(PAR22-AN2-AN2) D2(K) = (AN-0.1D+01)*(AN-0.2D+01)*PAR2 D1(K+1) = (AN+0.3D+01)*(AN+0.4D+01)*PAR2 V(K+2) = AC+(AN2-0.4D+01)*AS AN = AN+0.2D+01 70 CONTINUE AN2 = AN*AN D(NOEQU) = -0.2D+01*(AN2-0.4D+01)*(PAR22-AN2-AN2) V(NOEQU+2) = AC+(AN2-0.4D+01)*AS V(3) = V(3)-0.42D+02*PAR2*V(2) ASS = PARINT*COSPAR ASAP = (((((0.105D+03*PAR2-0.63D+02)*ASS+(0.210D+03*PAR2 1 -0.1D+01)*SINPAR)/AN2+(0.15D+02*PAR2-0.1D+01)*SINPAR- 2 0.15D+02*ASS)/AN2-0.3D+01*ASS-SINPAR)/AN2-SINPAR)/AN2 V(NOEQU+2) = V(NOEQU+2)-0.2D+01*ASAP*PAR2*(AN-0.1D+01) 1 *(AN-0.2D+01) C C SOLVE THE TRIDIAGONAL SYSTEM BY MEANS OF GAUSSIAN C ELIMINATION WITH PARTIAL PIVOTING. C C *** CALL TO DGTSL MUST BE REPLACED BY CALL TO C *** DOUBLE PRECISION VERSION OF LINPACK ROUTINE SGTSL C CALL DGTSL(NOEQU,D1,D,D2,V(3),IERS) GO TO 100 C C COMPUTE THE CHEBYSHEV MOMENTS BY MEANS OF FORWARD RECURSION. C 80 AN = 0.3D+01 DO 90 I = 3,12 AN2 = AN*AN V(I) = ((AN2-0.4D+01)*(0.2D+01*(PAR22-AN2-AN2)*V(I-1)+AS) 1 +AC-PAR2*(AN+0.1D+01)*(AN+0.2D+01)*V(I-2)) 2 /(PAR2*(AN-0.1D+01)*(AN-0.2D+01)) AN = AN+0.2D+01 90 CONTINUE 100 DO 110 J = 1,12 CHEBMO(M,2*J) = V(J) 110 CONTINUE 120 IF (NRMOM.LT.MOMCOM) M = NRMOM+1 IF (MOMCOM.LT.(MAXP1-1).AND.NRMOM.GE.MOMCOM) MOMCOM = MOMCOM+1 C C COMPUTE THE COEFFICIENTS OF THE CHEBYSHEV EXPANSIONS C OF DEGREES 12 AND 24 OF THE FUNCTION F. C FVAL(1) = 0.5D+00*F(CENTR+HLGTH) FVAL(13) = F(CENTR) FVAL(25) = 0.5D+00*F(CENTR-HLGTH) DO 130 I = 2,12 ISYM = 26-I FVAL(I) = F(HLGTH*X(I-1)+CENTR) FVAL(ISYM) = F(CENTR-HLGTH*X(I-1)) 130 CONTINUE CALL DQCHEB(X,FVAL,CHEB12,CHEB24) C C COMPUTE THE INTEGRAL AND ERROR ESTIMATES. C RESC12 = CHEB12(13)*CHEBMO(M,13) RESS12 = 0.0D+00 K = 11 DO 140 J = 1,6 RESC12 = RESC12+CHEB12(K)*CHEBMO(M,K) RESS12 = RESS12+CHEB12(K+1)*CHEBMO(M,K+1) K = K-2 140 CONTINUE RESC24 = CHEB24(25)*CHEBMO(M,25) RESS24 = 0.0D+00 RESABS = ABS(CHEB24(25)) K = 23 DO 150 J = 1,12 RESC24 = RESC24+CHEB24(K)*CHEBMO(M,K) RESS24 = RESS24+CHEB24(K+1)*CHEBMO(M,K+1) RESABS = ABS(CHEB24(K))+ABS(CHEB24(K+1)) K = K-2 150 CONTINUE ESTC = ABS(RESC24-RESC12) ESTS = ABS(RESS24-RESS12) RESABS = RESABS*ABS(HLGTH) IF(INTEGR.EQ.2) GO TO 160 RESULT = CONC*RESC24-CONS*RESS24 ABSERR = ABS(CONC*ESTC)+ABS(CONS*ESTS) GO TO 170 160 RESULT = CONC*RESS24+CONS*RESC24 ABSERR = ABS(CONC*ESTS)+ABS(CONS*ESTC) 170 RETURN END