1*DECK DQC25F 2 SUBROUTINE DQC25F (F, A, B, OMEGA, INTEGR, NRMOM, MAXP1, KSAVE, 3 + RESULT, ABSERR, NEVAL, RESABS, RESASC, MOMCOM, CHEBMO) 4C***BEGIN PROLOGUE DQC25F 5C***PURPOSE To compute the integral I=Integral of F(X) over (A,B) 6C Where W(X) = COS(OMEGA*X) or W(X)=SIN(OMEGA*X) and to 7C compute J = Integral of ABS(F) over (A,B). For small value 8C of OMEGA or small intervals (A,B) the 15-point GAUSS-KRONRO 9C Rule is used. Otherwise a generalized CLENSHAW-CURTIS 10C method is used. 11C***LIBRARY SLATEC (QUADPACK) 12C***CATEGORY H2A2A2 13C***TYPE DOUBLE PRECISION (QC25F-S, DQC25F-D) 14C***KEYWORDS CLENSHAW-CURTIS METHOD, GAUSS-KRONROD RULES, 15C INTEGRATION RULES FOR FUNCTIONS WITH COS OR SIN FACTOR, 16C QUADPACK, QUADRATURE 17C***AUTHOR Piessens, Robert 18C Applied Mathematics and Programming Division 19C K. U. Leuven 20C de Doncker, Elise 21C Applied Mathematics and Programming Division 22C K. U. Leuven 23C***DESCRIPTION 24C 25C Integration rules for functions with COS or SIN factor 26C Standard fortran subroutine 27C Double precision version 28C 29C PARAMETERS 30C ON ENTRY 31C F - Double precision 32C Function subprogram defining the integrand 33C function F(X). The actual name for F needs to 34C be declared E X T E R N A L in the calling program. 35C 36C A - Double precision 37C Lower limit of integration 38C 39C B - Double precision 40C Upper limit of integration 41C 42C OMEGA - Double precision 43C Parameter in the WEIGHT function 44C 45C INTEGR - Integer 46C Indicates which WEIGHT function is to be used 47C INTEGR = 1 W(X) = COS(OMEGA*X) 48C INTEGR = 2 W(X) = SIN(OMEGA*X) 49C 50C NRMOM - Integer 51C The length of interval (A,B) is equal to the length 52C of the original integration interval divided by 53C 2**NRMOM (we suppose that the routine is used in an 54C adaptive integration process, otherwise set 55C NRMOM = 0). NRMOM must be zero at the first call. 56C 57C MAXP1 - Integer 58C Gives an upper bound on the number of Chebyshev 59C moments which can be stored, i.e. for the 60C intervals of lengths ABS(BB-AA)*2**(-L), 61C L = 0,1,2, ..., MAXP1-2. 62C 63C KSAVE - Integer 64C Key which is one when the moments for the 65C current interval have been computed 66C 67C ON RETURN 68C RESULT - Double precision 69C Approximation to the integral I 70C 71C ABSERR - Double precision 72C Estimate of the modulus of the absolute 73C error, which should equal or exceed ABS(I-RESULT) 74C 75C NEVAL - Integer 76C Number of integrand evaluations 77C 78C RESABS - Double precision 79C Approximation to the integral J 80C 81C RESASC - Double precision 82C Approximation to the integral of ABS(F-I/(B-A)) 83C 84C ON ENTRY AND RETURN 85C MOMCOM - Integer 86C For each interval length we need to compute the 87C Chebyshev moments. MOMCOM counts the number of 88C intervals for which these moments have already been 89C computed. If NRMOM.LT.MOMCOM or KSAVE = 1, the 90C Chebyshev moments for the interval (A,B) have 91C already been computed and stored, otherwise we 92C compute them and we increase MOMCOM. 93C 94C CHEBMO - Double precision 95C Array of dimension at least (MAXP1,25) containing 96C the modified Chebyshev moments for the first MOMCOM 97C MOMCOM interval lengths 98C 99C ...................................................................... 100C 101C***REFERENCES (NONE) 102C***ROUTINES CALLED D1MACH, DGTSL, DQCHEB, DQK15W, DQWGTF 103C***REVISION HISTORY (YYMMDD) 104C 810101 DATE WRITTEN 105C 890531 Changed all specific intrinsics to generic. (WRB) 106C 890531 REVISION DATE from Version 3.2 107C 891214 Prologue converted to Version 4.0 format. (BAB) 108C***END PROLOGUE DQC25F 109C 110 DOUBLE PRECISION A,ABSERR,AC,AN,AN2,AS,ASAP,ASS,B,CENTR,CHEBMO, 111 1 CHEB12,CHEB24,CONC,CONS,COSPAR,D,DQWGTF,D1, 112 2 D1MACH,D2,ESTC,ESTS,F,FVAL,HLGTH,OFLOW,OMEGA,PARINT,PAR2,PAR22, 113 3 P2,P3,P4,RESABS,RESASC,RESC12,RESC24,RESS12,RESS24,RESULT, 114 4 SINPAR,V,X 115 INTEGER I,IERS,INTEGR,ISYM,J,K,KSAVE,M,MOMCOM,NEVAL,MAXP1, 116 1 NOEQU,NOEQ1,NRMOM 117C 118 DIMENSION CHEBMO(MAXP1,25),CHEB12(13),CHEB24(25),D(25),D1(25), 119 1 D2(25),FVAL(25),V(28),X(11) 120C 121 EXTERNAL F, DQWGTF 122C 123C THE VECTOR X CONTAINS THE VALUES COS(K*PI/24) 124C K = 1, ...,11, TO BE USED FOR THE CHEBYSHEV EXPANSION OF F 125C 126 SAVE X 127 DATA X(1) / 0.9914448613 7381041114 4557526928 563D0 / 128 DATA X(2) / 0.9659258262 8906828674 9743199728 897D0 / 129 DATA X(3) / 0.9238795325 1128675612 8183189396 788D0 / 130 DATA X(4) / 0.8660254037 8443864676 3723170752 936D0 / 131 DATA X(5) / 0.7933533402 9123516457 9776961501 299D0 / 132 DATA X(6) / 0.7071067811 8654752440 0844362104 849D0 / 133 DATA X(7) / 0.6087614290 0872063941 6097542898 164D0 / 134 DATA X(8) / 0.5000000000 0000000000 0000000000 000D0 / 135 DATA X(9) / 0.3826834323 6508977172 8459984030 399D0 / 136 DATA X(10) / 0.2588190451 0252076234 8898837624 048D0 / 137 DATA X(11) / 0.1305261922 2005159154 8406227895 489D0 / 138C 139C LIST OF MAJOR VARIABLES 140C ----------------------- 141C 142C CENTR - MID POINT OF THE INTEGRATION INTERVAL 143C HLGTH - HALF-LENGTH OF THE INTEGRATION INTERVAL 144C FVAL - VALUE OF THE FUNCTION F AT THE POINTS 145C (B-A)*0.5*COS(K*PI/12) + (B+A)*0.5, K = 0, ..., 24 146C CHEB12 - COEFFICIENTS OF THE CHEBYSHEV SERIES EXPANSION 147C OF DEGREE 12, FOR THE FUNCTION F, IN THE 148C INTERVAL (A,B) 149C CHEB24 - COEFFICIENTS OF THE CHEBYSHEV SERIES EXPANSION 150C OF DEGREE 24, FOR THE FUNCTION F, IN THE 151C INTERVAL (A,B) 152C RESC12 - APPROXIMATION TO THE INTEGRAL OF 153C COS(0.5*(B-A)*OMEGA*X)*F(0.5*(B-A)*X+0.5*(B+A)) 154C OVER (-1,+1), USING THE CHEBYSHEV SERIES 155C EXPANSION OF DEGREE 12 156C RESC24 - APPROXIMATION TO THE SAME INTEGRAL, USING THE 157C CHEBYSHEV SERIES EXPANSION OF DEGREE 24 158C RESS12 - THE ANALOGUE OF RESC12 FOR THE SINE 159C RESS24 - THE ANALOGUE OF RESC24 FOR THE SINE 160C 161C 162C MACHINE DEPENDENT CONSTANT 163C -------------------------- 164C 165C OFLOW IS THE LARGEST POSITIVE MAGNITUDE. 166C 167C***FIRST EXECUTABLE STATEMENT DQC25F 168 OFLOW = D1MACH(2) 169C 170 CENTR = 0.5D+00*(B+A) 171 HLGTH = 0.5D+00*(B-A) 172 PARINT = OMEGA*HLGTH 173C 174C COMPUTE THE INTEGRAL USING THE 15-POINT GAUSS-KRONROD 175C FORMULA IF THE VALUE OF THE PARAMETER IN THE INTEGRAND 176C IS SMALL. 177C 178 IF(ABS(PARINT).GT.0.2D+01) GO TO 10 179 CALL DQK15W(F,DQWGTF,OMEGA,P2,P3,P4,INTEGR,A,B,RESULT, 180 1 ABSERR,RESABS,RESASC) 181 NEVAL = 15 182 GO TO 170 183C 184C COMPUTE THE INTEGRAL USING THE GENERALIZED CLENSHAW- 185C CURTIS METHOD. 186C 187 10 CONC = HLGTH*COS(CENTR*OMEGA) 188 CONS = HLGTH*SIN(CENTR*OMEGA) 189 RESASC = OFLOW 190 NEVAL = 25 191C 192C CHECK WHETHER THE CHEBYSHEV MOMENTS FOR THIS INTERVAL 193C HAVE ALREADY BEEN COMPUTED. 194C 195 IF(NRMOM.LT.MOMCOM.OR.KSAVE.EQ.1) GO TO 120 196C 197C COMPUTE A NEW SET OF CHEBYSHEV MOMENTS. 198C 199 M = MOMCOM+1 200 PAR2 = PARINT*PARINT 201 PAR22 = PAR2+0.2D+01 202 SINPAR = SIN(PARINT) 203 COSPAR = COS(PARINT) 204C 205C COMPUTE THE CHEBYSHEV MOMENTS WITH RESPECT TO COSINE. 206C 207 V(1) = 0.2D+01*SINPAR/PARINT 208 V(2) = (0.8D+01*COSPAR+(PAR2+PAR2-0.8D+01)*SINPAR/PARINT)/PAR2 209 V(3) = (0.32D+02*(PAR2-0.12D+02)*COSPAR+(0.2D+01* 210 1 ((PAR2-0.80D+02)*PAR2+0.192D+03)*SINPAR)/PARINT)/(PAR2*PAR2) 211 AC = 0.8D+01*COSPAR 212 AS = 0.24D+02*PARINT*SINPAR 213 IF(ABS(PARINT).GT.0.24D+02) GO TO 30 214C 215C COMPUTE THE CHEBYSHEV MOMENTS AS THE SOLUTIONS OF A 216C BOUNDARY VALUE PROBLEM WITH 1 INITIAL VALUE (V(3)) AND 1 217C END VALUE (COMPUTED USING AN ASYMPTOTIC FORMULA). 218C 219 NOEQU = 25 220 NOEQ1 = NOEQU-1 221 AN = 0.6D+01 222 DO 20 K = 1,NOEQ1 223 AN2 = AN*AN 224 D(K) = -0.2D+01*(AN2-0.4D+01)*(PAR22-AN2-AN2) 225 D2(K) = (AN-0.1D+01)*(AN-0.2D+01)*PAR2 226 D1(K+1) = (AN+0.3D+01)*(AN+0.4D+01)*PAR2 227 V(K+3) = AS-(AN2-0.4D+01)*AC 228 AN = AN+0.2D+01 229 20 CONTINUE 230 AN2 = AN*AN 231 D(NOEQU) = -0.2D+01*(AN2-0.4D+01)*(PAR22-AN2-AN2) 232 V(NOEQU+3) = AS-(AN2-0.4D+01)*AC 233 V(4) = V(4)-0.56D+02*PAR2*V(3) 234 ASS = PARINT*SINPAR 235 ASAP = (((((0.210D+03*PAR2-0.1D+01)*COSPAR-(0.105D+03*PAR2 236 1 -0.63D+02)*ASS)/AN2-(0.1D+01-0.15D+02*PAR2)*COSPAR 237 2 +0.15D+02*ASS)/AN2-COSPAR+0.3D+01*ASS)/AN2-COSPAR)/AN2 238 V(NOEQU+3) = V(NOEQU+3)-0.2D+01*ASAP*PAR2*(AN-0.1D+01)* 239 1 (AN-0.2D+01) 240C 241C SOLVE THE TRIDIAGONAL SYSTEM BY MEANS OF GAUSSIAN 242C ELIMINATION WITH PARTIAL PIVOTING. 243C 244C *** CALL TO DGTSL MUST BE REPLACED BY CALL TO 245C *** DOUBLE PRECISION VERSION OF LINPACK ROUTINE SGTSL 246C 247 CALL DGTSL(NOEQU,D1,D,D2,V(4),IERS) 248 GO TO 50 249C 250C COMPUTE THE CHEBYSHEV MOMENTS BY MEANS OF FORWARD 251C RECURSION. 252C 253 30 AN = 0.4D+01 254 DO 40 I = 4,13 255 AN2 = AN*AN 256 V(I) = ((AN2-0.4D+01)*(0.2D+01*(PAR22-AN2-AN2)*V(I-1)-AC) 257 1 +AS-PAR2*(AN+0.1D+01)*(AN+0.2D+01)*V(I-2))/ 258 2 (PAR2*(AN-0.1D+01)*(AN-0.2D+01)) 259 AN = AN+0.2D+01 260 40 CONTINUE 261 50 DO 60 J = 1,13 262 CHEBMO(M,2*J-1) = V(J) 263 60 CONTINUE 264C 265C COMPUTE THE CHEBYSHEV MOMENTS WITH RESPECT TO SINE. 266C 267 V(1) = 0.2D+01*(SINPAR-PARINT*COSPAR)/PAR2 268 V(2) = (0.18D+02-0.48D+02/PAR2)*SINPAR/PAR2 269 1 +(-0.2D+01+0.48D+02/PAR2)*COSPAR/PARINT 270 AC = -0.24D+02*PARINT*COSPAR 271 AS = -0.8D+01*SINPAR 272 IF(ABS(PARINT).GT.0.24D+02) GO TO 80 273C 274C COMPUTE THE CHEBYSHEV MOMENTS AS THE SOLUTIONS OF A BOUNDARY 275C VALUE PROBLEM WITH 1 INITIAL VALUE (V(2)) AND 1 END VALUE 276C (COMPUTED USING AN ASYMPTOTIC FORMULA). 277C 278 AN = 0.5D+01 279 DO 70 K = 1,NOEQ1 280 AN2 = AN*AN 281 D(K) = -0.2D+01*(AN2-0.4D+01)*(PAR22-AN2-AN2) 282 D2(K) = (AN-0.1D+01)*(AN-0.2D+01)*PAR2 283 D1(K+1) = (AN+0.3D+01)*(AN+0.4D+01)*PAR2 284 V(K+2) = AC+(AN2-0.4D+01)*AS 285 AN = AN+0.2D+01 286 70 CONTINUE 287 AN2 = AN*AN 288 D(NOEQU) = -0.2D+01*(AN2-0.4D+01)*(PAR22-AN2-AN2) 289 V(NOEQU+2) = AC+(AN2-0.4D+01)*AS 290 V(3) = V(3)-0.42D+02*PAR2*V(2) 291 ASS = PARINT*COSPAR 292 ASAP = (((((0.105D+03*PAR2-0.63D+02)*ASS+(0.210D+03*PAR2 293 1 -0.1D+01)*SINPAR)/AN2+(0.15D+02*PAR2-0.1D+01)*SINPAR- 294 2 0.15D+02*ASS)/AN2-0.3D+01*ASS-SINPAR)/AN2-SINPAR)/AN2 295 V(NOEQU+2) = V(NOEQU+2)-0.2D+01*ASAP*PAR2*(AN-0.1D+01) 296 1 *(AN-0.2D+01) 297C 298C SOLVE THE TRIDIAGONAL SYSTEM BY MEANS OF GAUSSIAN 299C ELIMINATION WITH PARTIAL PIVOTING. 300C 301C *** CALL TO DGTSL MUST BE REPLACED BY CALL TO 302C *** DOUBLE PRECISION VERSION OF LINPACK ROUTINE SGTSL 303C 304 CALL DGTSL(NOEQU,D1,D,D2,V(3),IERS) 305 GO TO 100 306C 307C COMPUTE THE CHEBYSHEV MOMENTS BY MEANS OF FORWARD RECURSION. 308C 309 80 AN = 0.3D+01 310 DO 90 I = 3,12 311 AN2 = AN*AN 312 V(I) = ((AN2-0.4D+01)*(0.2D+01*(PAR22-AN2-AN2)*V(I-1)+AS) 313 1 +AC-PAR2*(AN+0.1D+01)*(AN+0.2D+01)*V(I-2)) 314 2 /(PAR2*(AN-0.1D+01)*(AN-0.2D+01)) 315 AN = AN+0.2D+01 316 90 CONTINUE 317 100 DO 110 J = 1,12 318 CHEBMO(M,2*J) = V(J) 319 110 CONTINUE 320 120 IF (NRMOM.LT.MOMCOM) M = NRMOM+1 321 IF (MOMCOM.LT.(MAXP1-1).AND.NRMOM.GE.MOMCOM) MOMCOM = MOMCOM+1 322C 323C COMPUTE THE COEFFICIENTS OF THE CHEBYSHEV EXPANSIONS 324C OF DEGREES 12 AND 24 OF THE FUNCTION F. 325C 326 FVAL(1) = 0.5D+00*F(CENTR+HLGTH) 327 FVAL(13) = F(CENTR) 328 FVAL(25) = 0.5D+00*F(CENTR-HLGTH) 329 DO 130 I = 2,12 330 ISYM = 26-I 331 FVAL(I) = F(HLGTH*X(I-1)+CENTR) 332 FVAL(ISYM) = F(CENTR-HLGTH*X(I-1)) 333 130 CONTINUE 334 CALL DQCHEB(X,FVAL,CHEB12,CHEB24) 335C 336C COMPUTE THE INTEGRAL AND ERROR ESTIMATES. 337C 338 RESC12 = CHEB12(13)*CHEBMO(M,13) 339 RESS12 = 0.0D+00 340 K = 11 341 DO 140 J = 1,6 342 RESC12 = RESC12+CHEB12(K)*CHEBMO(M,K) 343 RESS12 = RESS12+CHEB12(K+1)*CHEBMO(M,K+1) 344 K = K-2 345 140 CONTINUE 346 RESC24 = CHEB24(25)*CHEBMO(M,25) 347 RESS24 = 0.0D+00 348 RESABS = ABS(CHEB24(25)) 349 K = 23 350 DO 150 J = 1,12 351 RESC24 = RESC24+CHEB24(K)*CHEBMO(M,K) 352 RESS24 = RESS24+CHEB24(K+1)*CHEBMO(M,K+1) 353 RESABS = RESABS + ABS(CHEB24(K))+ABS(CHEB24(K+1)) 354 K = K-2 355 150 CONTINUE 356 ESTC = ABS(RESC24-RESC12) 357 ESTS = ABS(RESS24-RESS12) 358 RESABS = RESABS*ABS(HLGTH) 359 IF(INTEGR.EQ.2) GO TO 160 360 RESULT = CONC*RESC24-CONS*RESS24 361 ABSERR = ABS(CONC*ESTC)+ABS(CONS*ESTS) 362 GO TO 170 363 160 RESULT = CONC*RESS24+CONS*RESC24 364 ABSERR = ABS(CONC*ESTS)+ABS(CONS*ESTC) 365 170 RETURN 366 END 367