2
SUBROUTINE DQC25F (F, A, B, OMEGA, INTEGR, NRMOM, MAXP1, KSAVE,
3
+ RESULT, ABSERR, NEVAL, RESABS, RESASC, MOMCOM, CHEBMO)
4
C***BEGIN PROLOGUE DQC25F
5
C***PURPOSE To compute the integral I=Integral of F(X) over (A,B)
6
C Where W(X) = COS(OMEGA*X) or W(X)=SIN(OMEGA*X) and to
7
C compute J = Integral of ABS(F) over (A,B). For small value
8
C of OMEGA or small intervals (A,B) the 15-point GAUSS-KRONRO
9
C Rule is used. Otherwise a generalized CLENSHAW-CURTIS
11
C***LIBRARY SLATEC (QUADPACK)
13
C***TYPE DOUBLE PRECISION (QC25F-S, DQC25F-D)
14
C***KEYWORDS CLENSHAW-CURTIS METHOD, GAUSS-KRONROD RULES,
15
C INTEGRATION RULES FOR FUNCTIONS WITH COS OR SIN FACTOR,
16
C QUADPACK, QUADRATURE
17
C***AUTHOR Piessens, Robert
18
C Applied Mathematics and Programming Division
21
C Applied Mathematics and Programming Division
25
C Integration rules for functions with COS or SIN factor
26
C Standard fortran subroutine
27
C Double precision version
31
C F - Double precision
32
C Function subprogram defining the integrand
33
C function F(X). The actual name for F needs to
34
C be declared E X T E R N A L in the calling program.
36
C A - Double precision
37
C Lower limit of integration
39
C B - Double precision
40
C Upper limit of integration
42
C OMEGA - Double precision
43
C Parameter in the WEIGHT function
46
C Indicates which WEIGHT function is to be used
47
C INTEGR = 1 W(X) = COS(OMEGA*X)
48
C INTEGR = 2 W(X) = SIN(OMEGA*X)
51
C The length of interval (A,B) is equal to the length
52
C of the original integration interval divided by
53
C 2**NRMOM (we suppose that the routine is used in an
54
C adaptive integration process, otherwise set
55
C NRMOM = 0). NRMOM must be zero at the first call.
58
C Gives an upper bound on the number of Chebyshev
59
C moments which can be stored, i.e. for the
60
C intervals of lengths ABS(BB-AA)*2**(-L),
61
C L = 0,1,2, ..., MAXP1-2.
64
C Key which is one when the moments for the
65
C current interval have been computed
68
C RESULT - Double precision
69
C Approximation to the integral I
71
C ABSERR - Double precision
72
C Estimate of the modulus of the absolute
73
C error, which should equal or exceed ABS(I-RESULT)
76
C Number of integrand evaluations
78
C RESABS - Double precision
79
C Approximation to the integral J
81
C RESASC - Double precision
82
C Approximation to the integral of ABS(F-I/(B-A))
86
C For each interval length we need to compute the
87
C Chebyshev moments. MOMCOM counts the number of
88
C intervals for which these moments have already been
89
C computed. If NRMOM.LT.MOMCOM or KSAVE = 1, the
90
C Chebyshev moments for the interval (A,B) have
91
C already been computed and stored, otherwise we
92
C compute them and we increase MOMCOM.
94
C CHEBMO - Double precision
95
C Array of dimension at least (MAXP1,25) containing
96
C the modified Chebyshev moments for the first MOMCOM
97
C MOMCOM interval lengths
99
C ......................................................................
101
C***REFERENCES (NONE)
102
C***ROUTINES CALLED D1MACH, DGTSL, DQCHEB, DQK15W, DQWGTF
103
C***REVISION HISTORY (YYMMDD)
104
C 810101 DATE WRITTEN
105
C 890531 Changed all specific intrinsics to generic. (WRB)
106
C 890531 REVISION DATE from Version 3.2
107
C 891214 Prologue converted to Version 4.0 format. (BAB)
108
C***END PROLOGUE DQC25F
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,
115
INTEGER I,IERS,INTEGR,ISYM,J,K,KSAVE,M,MOMCOM,NEVAL,MAXP1,
118
DIMENSION CHEBMO(MAXP1,25),CHEB12(13),CHEB24(25),D(25),D1(25),
119
1 D2(25),FVAL(25),V(28),X(11)
123
C THE VECTOR X CONTAINS THE VALUES COS(K*PI/24)
124
C K = 1, ...,11, TO BE USED FOR THE CHEBYSHEV EXPANSION OF F
127
DATA X(1),X(2),X(3),X(4),X(5),X(6),X(7),X(8),X(9),X(10),X(11)/
128
1 0.9914448613738104D+00, 0.9659258262890683D+00,
129
2 0.9238795325112868D+00, 0.8660254037844386D+00,
130
3 0.7933533402912352D+00, 0.7071067811865475D+00,
131
4 0.6087614290087206D+00, 0.5000000000000000D+00,
132
5 0.3826834323650898D+00, 0.2588190451025208D+00,
133
6 0.1305261922200516D+00/
135
C LIST OF MAJOR VARIABLES
136
C -----------------------
138
C CENTR - MID POINT OF THE INTEGRATION INTERVAL
139
C HLGTH - HALF-LENGTH OF THE INTEGRATION INTERVAL
140
C FVAL - VALUE OF THE FUNCTION F AT THE POINTS
141
C (B-A)*0.5*COS(K*PI/12) + (B+A)*0.5, K = 0, ..., 24
142
C CHEB12 - COEFFICIENTS OF THE CHEBYSHEV SERIES EXPANSION
143
C OF DEGREE 12, FOR THE FUNCTION F, IN THE
145
C CHEB24 - COEFFICIENTS OF THE CHEBYSHEV SERIES EXPANSION
146
C OF DEGREE 24, FOR THE FUNCTION F, IN THE
148
C RESC12 - APPROXIMATION TO THE INTEGRAL OF
149
C COS(0.5*(B-A)*OMEGA*X)*F(0.5*(B-A)*X+0.5*(B+A))
150
C OVER (-1,+1), USING THE CHEBYSHEV SERIES
151
C EXPANSION OF DEGREE 12
152
C RESC24 - APPROXIMATION TO THE SAME INTEGRAL, USING THE
153
C CHEBYSHEV SERIES EXPANSION OF DEGREE 24
154
C RESS12 - THE ANALOGUE OF RESC12 FOR THE SINE
155
C RESS24 - THE ANALOGUE OF RESC24 FOR THE SINE
158
C MACHINE DEPENDENT CONSTANT
159
C --------------------------
161
C OFLOW IS THE LARGEST POSITIVE MAGNITUDE.
163
C***FIRST EXECUTABLE STATEMENT DQC25F
166
CENTR = 0.5D+00*(B+A)
167
HLGTH = 0.5D+00*(B-A)
170
C COMPUTE THE INTEGRAL USING THE 15-POINT GAUSS-KRONROD
171
C FORMULA IF THE VALUE OF THE PARAMETER IN THE INTEGRAND
174
IF(ABS(PARINT).GT.0.2D+01) GO TO 10
175
CALL DQK15W(F,DQWGTF,OMEGA,P2,P3,P4,INTEGR,A,B,RESULT,
176
1 ABSERR,RESABS,RESASC)
180
C COMPUTE THE INTEGRAL USING THE GENERALIZED CLENSHAW-
183
10 CONC = HLGTH*COS(CENTR*OMEGA)
184
CONS = HLGTH*SIN(CENTR*OMEGA)
188
C CHECK WHETHER THE CHEBYSHEV MOMENTS FOR THIS INTERVAL
189
C HAVE ALREADY BEEN COMPUTED.
191
IF(NRMOM.LT.MOMCOM.OR.KSAVE.EQ.1) GO TO 120
193
C COMPUTE A NEW SET OF CHEBYSHEV MOMENTS.
201
C COMPUTE THE CHEBYSHEV MOMENTS WITH RESPECT TO COSINE.
203
V(1) = 0.2D+01*SINPAR/PARINT
204
V(2) = (0.8D+01*COSPAR+(PAR2+PAR2-0.8D+01)*SINPAR/PARINT)/PAR2
205
V(3) = (0.32D+02*(PAR2-0.12D+02)*COSPAR+(0.2D+01*
206
1 ((PAR2-0.80D+02)*PAR2+0.192D+03)*SINPAR)/PARINT)/(PAR2*PAR2)
208
AS = 0.24D+02*PARINT*SINPAR
209
IF(ABS(PARINT).GT.0.24D+02) GO TO 30
211
C COMPUTE THE CHEBYSHEV MOMENTS AS THE SOLUTIONS OF A
212
C BOUNDARY VALUE PROBLEM WITH 1 INITIAL VALUE (V(3)) AND 1
213
C END VALUE (COMPUTED USING AN ASYMPTOTIC FORMULA).
220
D(K) = -0.2D+01*(AN2-0.4D+01)*(PAR22-AN2-AN2)
221
D2(K) = (AN-0.1D+01)*(AN-0.2D+01)*PAR2
222
D1(K+1) = (AN+0.3D+01)*(AN+0.4D+01)*PAR2
223
V(K+3) = AS-(AN2-0.4D+01)*AC
227
D(NOEQU) = -0.2D+01*(AN2-0.4D+01)*(PAR22-AN2-AN2)
228
V(NOEQU+3) = AS-(AN2-0.4D+01)*AC
229
V(4) = V(4)-0.56D+02*PAR2*V(3)
231
ASAP = (((((0.210D+03*PAR2-0.1D+01)*COSPAR-(0.105D+03*PAR2
232
1 -0.63D+02)*ASS)/AN2-(0.1D+01-0.15D+02*PAR2)*COSPAR
233
2 +0.15D+02*ASS)/AN2-COSPAR+0.3D+01*ASS)/AN2-COSPAR)/AN2
234
V(NOEQU+3) = V(NOEQU+3)-0.2D+01*ASAP*PAR2*(AN-0.1D+01)*
237
C SOLVE THE TRIDIAGONAL SYSTEM BY MEANS OF GAUSSIAN
238
C ELIMINATION WITH PARTIAL PIVOTING.
240
C *** CALL TO DGTSL MUST BE REPLACED BY CALL TO
241
C *** DOUBLE PRECISION VERSION OF LINPACK ROUTINE SGTSL
243
CALL DGTSL(NOEQU,D1,D,D2,V(4),IERS)
246
C COMPUTE THE CHEBYSHEV MOMENTS BY MEANS OF FORWARD
252
V(I) = ((AN2-0.4D+01)*(0.2D+01*(PAR22-AN2-AN2)*V(I-1)-AC)
253
1 +AS-PAR2*(AN+0.1D+01)*(AN+0.2D+01)*V(I-2))/
254
2 (PAR2*(AN-0.1D+01)*(AN-0.2D+01))
258
CHEBMO(M,2*J-1) = V(J)
261
C COMPUTE THE CHEBYSHEV MOMENTS WITH RESPECT TO SINE.
263
V(1) = 0.2D+01*(SINPAR-PARINT*COSPAR)/PAR2
264
V(2) = (0.18D+02-0.48D+02/PAR2)*SINPAR/PAR2
265
1 +(-0.2D+01+0.48D+02/PAR2)*COSPAR/PARINT
266
AC = -0.24D+02*PARINT*COSPAR
268
IF(ABS(PARINT).GT.0.24D+02) GO TO 80
270
C COMPUTE THE CHEBYSHEV MOMENTS AS THE SOLUTIONS OF A BOUNDARY
271
C VALUE PROBLEM WITH 1 INITIAL VALUE (V(2)) AND 1 END VALUE
272
C (COMPUTED USING AN ASYMPTOTIC FORMULA).
277
D(K) = -0.2D+01*(AN2-0.4D+01)*(PAR22-AN2-AN2)
278
D2(K) = (AN-0.1D+01)*(AN-0.2D+01)*PAR2
279
D1(K+1) = (AN+0.3D+01)*(AN+0.4D+01)*PAR2
280
V(K+2) = AC+(AN2-0.4D+01)*AS
284
D(NOEQU) = -0.2D+01*(AN2-0.4D+01)*(PAR22-AN2-AN2)
285
V(NOEQU+2) = AC+(AN2-0.4D+01)*AS
286
V(3) = V(3)-0.42D+02*PAR2*V(2)
288
ASAP = (((((0.105D+03*PAR2-0.63D+02)*ASS+(0.210D+03*PAR2
289
1 -0.1D+01)*SINPAR)/AN2+(0.15D+02*PAR2-0.1D+01)*SINPAR-
290
2 0.15D+02*ASS)/AN2-0.3D+01*ASS-SINPAR)/AN2-SINPAR)/AN2
291
V(NOEQU+2) = V(NOEQU+2)-0.2D+01*ASAP*PAR2*(AN-0.1D+01)
294
C SOLVE THE TRIDIAGONAL SYSTEM BY MEANS OF GAUSSIAN
295
C ELIMINATION WITH PARTIAL PIVOTING.
297
C *** CALL TO DGTSL MUST BE REPLACED BY CALL TO
298
C *** DOUBLE PRECISION VERSION OF LINPACK ROUTINE SGTSL
300
CALL DGTSL(NOEQU,D1,D,D2,V(3),IERS)
303
C COMPUTE THE CHEBYSHEV MOMENTS BY MEANS OF FORWARD RECURSION.
308
V(I) = ((AN2-0.4D+01)*(0.2D+01*(PAR22-AN2-AN2)*V(I-1)+AS)
309
1 +AC-PAR2*(AN+0.1D+01)*(AN+0.2D+01)*V(I-2))
310
2 /(PAR2*(AN-0.1D+01)*(AN-0.2D+01))
316
120 IF (NRMOM.LT.MOMCOM) M = NRMOM+1
317
IF (MOMCOM.LT.(MAXP1-1).AND.NRMOM.GE.MOMCOM) MOMCOM = MOMCOM+1
319
C COMPUTE THE COEFFICIENTS OF THE CHEBYSHEV EXPANSIONS
320
C OF DEGREES 12 AND 24 OF THE FUNCTION F.
322
FVAL(1) = 0.5D+00*F(CENTR+HLGTH)
324
FVAL(25) = 0.5D+00*F(CENTR-HLGTH)
327
FVAL(I) = F(HLGTH*X(I-1)+CENTR)
328
FVAL(ISYM) = F(CENTR-HLGTH*X(I-1))
330
CALL DQCHEB(X,FVAL,CHEB12,CHEB24)
332
C COMPUTE THE INTEGRAL AND ERROR ESTIMATES.
334
RESC12 = CHEB12(13)*CHEBMO(M,13)
338
RESC12 = RESC12+CHEB12(K)*CHEBMO(M,K)
339
RESS12 = RESS12+CHEB12(K+1)*CHEBMO(M,K+1)
342
RESC24 = CHEB24(25)*CHEBMO(M,25)
344
RESABS = ABS(CHEB24(25))
347
RESC24 = RESC24+CHEB24(K)*CHEBMO(M,K)
348
RESS24 = RESS24+CHEB24(K+1)*CHEBMO(M,K+1)
349
RESABS = ABS(CHEB24(K))+ABS(CHEB24(K+1))
352
ESTC = ABS(RESC24-RESC12)
353
ESTS = ABS(RESS24-RESS12)
354
RESABS = RESABS*ABS(HLGTH)
355
IF(INTEGR.EQ.2) GO TO 160
356
RESULT = CONC*RESC24-CONS*RESS24
357
ABSERR = ABS(CONC*ESTC)+ABS(CONS*ESTS)
359
160 RESULT = CONC*RESS24+CONS*RESC24
360
ABSERR = ABS(CONC*ESTS)+ABS(CONS*ESTC)