|
|
SUBROUTINE QC25F(F,A,B,OMEGA,INTEGR,NRMOM,MAXP1,KSAVE,RESULT,
|
|
|
* ABSERR,NEVAL,RESABS,RESASC,MOMCOM,CHEBMO)
|
|
|
C***BEGIN PROLOGUE QC25F
|
|
|
C***DATE WRITTEN 810101 (YYMMDD)
|
|
|
C***REVISION DATE 830518 (YYMMDD)
|
|
|
C***CATEGORY NO. H2A2A2
|
|
|
C***KEYWORDS INTEGRATION RULES FOR FUNCTIONS WITH COS OR SIN
|
|
|
C FACTOR, CLENSHAW-CURTIS, GAUSS-KRONROD
|
|
|
C***AUTHOR PIESSENS,ROBERT,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN
|
|
|
C DE DONCKER,ELISE,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN
|
|
|
C***PURPOSE TO COMPUTE THE INTEGRAL I=INTEGRAL OF F(X) OVER (A,B)
|
|
|
C WHERE W(X) = COS(OMEGA*X) OR (WX)=SIN(OMEGA*X)
|
|
|
C AND TO COMPUTE J=INTEGRAL OF ABS(F) OVER (A,B). FOR SMALL
|
|
|
C VALUE OF OMEGA OR SMALL INTERVALS (A,B) 15-POINT GAUSS-
|
|
|
C KRONROD RULE USED. OTHERWISE GENERALIZED CLENSHAW-CURTIS US
|
|
|
C***DESCRIPTION
|
|
|
C
|
|
|
C INTEGRATION RULES FOR FUNCTIONS WITH COS OR SIN FACTOR
|
|
|
C STANDARD FORTRAN SUBROUTINE
|
|
|
C REAL VERSION
|
|
|
C
|
|
|
C PARAMETERS
|
|
|
C ON ENTRY
|
|
|
C F - REAL
|
|
|
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 - REAL
|
|
|
C LOWER LIMIT OF INTEGRATION
|
|
|
C
|
|
|
C B - REAL
|
|
|
C UPPER LIMIT OF INTEGRATION
|
|
|
C
|
|
|
C OMEGA - REAL
|
|
|
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 - REAL
|
|
|
C APPROXIMATION TO THE INTEGRAL I
|
|
|
C
|
|
|
C ABSERR - REAL
|
|
|
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 - REAL
|
|
|
C APPROXIMATION TO THE INTEGRAL J
|
|
|
C
|
|
|
C RESASC - REAL
|
|
|
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 - REAL
|
|
|
C ARRAY OF DIMENSION AT LEAST (MAXP1,25) CONTAINING
|
|
|
C THE MODIFIED CHEBYSHEV MOMENTS FOR THE FIRST MOMCOM
|
|
|
C MOMCOM INTERVAL LENGTHS
|
|
|
C
|
|
|
C***REFERENCES (NONE)
|
|
|
C***ROUTINES CALLED QCHEB,QK15W,QWGTF,R1MACH,SGTSL
|
|
|
C***END PROLOGUE QC25F
|
|
|
C
|
|
|
REAL A,ABSERR,AC,AN,AN2,AS,ASAP,ASS,B,CENTR,CHEBMO,
|
|
|
* CHEB12,CHEB24,CONC,CONS,COSPAR,D,QWGTF,
|
|
|
* D1,R1MACH,D2,ESTC,ESTS,F,FVAL,HLGTH,OFLOW,OMEGA,PARINT,PAR2,
|
|
|
* PAR22,P2,P3,P4,RESABS,RESASC,RESC12,RESC24,RESS12,RESS24,
|
|
|
* RESULT,SINPAR,V,X
|
|
|
INTEGER I,IERS,INTEGR,ISYM,J,K,KSAVE,M,MAXP1,MOMCOM,NEVAL,
|
|
|
* NOEQU,NOEQ1,NRMOM
|
|
|
C
|
|
|
DIMENSION CHEBMO(MAXP1,25),CHEB12(13),CHEB24(25),D(25),D1(25),
|
|
|
* D2(25),FVAL(25),V(28),X(11)
|
|
|
C
|
|
|
EXTERNAL F,QWGTF
|
|
|
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
|
|
|
DATA X(1),X(2),X(3),X(4),X(5),X(6),X(7),X(8),X(9),
|
|
|
* X(10),X(11)/
|
|
|
* 0.9914448613738104E+00, 0.9659258262890683E+00,
|
|
|
* 0.9238795325112868E+00, 0.8660254037844386E+00,
|
|
|
* 0.7933533402912352E+00, 0.7071067811865475E+00,
|
|
|
* 0.6087614290087206E+00, 0.5000000000000000E+00,
|
|
|
* 0.3826834323650898E+00, 0.2588190451025208E+00,
|
|
|
* 0.1305261922200516E+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,
|
|
|
C 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 QC25F
|
|
|
OFLOW = R1MACH(2)
|
|
|
C
|
|
|
CENTR = 0.5E+00*(B+A)
|
|
|
HLGTH = 0.5E+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.2E+01) GO TO 10
|
|
|
CALL QK15W(F,QWGTF,OMEGA,P2,P3,P4,INTEGR,A,B,RESULT,
|
|
|
* 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.2E+01
|
|
|
SINPAR = SIN(PARINT)
|
|
|
COSPAR = COS(PARINT)
|
|
|
C
|
|
|
C COMPUTE THE CHEBYSHEV MOMENTS WITH RESPECT TO COSINE.
|
|
|
C
|
|
|
V(1) = 0.2E+01*SINPAR/PARINT
|
|
|
V(2) = (0.8E+01*COSPAR+(PAR2+PAR2-0.8E+01)*SINPAR/
|
|
|
* PARINT)/PAR2
|
|
|
V(3) = (0.32E+02*(PAR2-0.12E+02)*COSPAR+(0.2E+01*
|
|
|
* ((PAR2-0.80E+02)*PAR2+0.192E+03)*SINPAR)/
|
|
|
* PARINT)/(PAR2*PAR2)
|
|
|
AC = 0.8E+01*COSPAR
|
|
|
AS = 0.24E+02*PARINT*SINPAR
|
|
|
IF(ABS(PARINT).GT.0.24E+02) GO TO 30
|
|
|
C
|
|
|
C COMPUTE THE CHEBYSHEV MOMENTS AS THE
|
|
|
C SOLUTIONS OF A BOUNDARY VALUE PROBLEM WITH 1
|
|
|
C INITIAL VALUE (V(3)) AND 1 END VALUE (COMPUTED
|
|
|
C USING AN ASYMPTOTIC FORMULA).
|
|
|
C
|
|
|
NOEQU = 25
|
|
|
NOEQ1 = NOEQU-1
|
|
|
AN = 0.6E+01
|
|
|
DO 20 K = 1,NOEQ1
|
|
|
AN2 = AN*AN
|
|
|
D(K) = -0.2E+01*(AN2-0.4E+01)*(PAR22-AN2-AN2)
|
|
|
D2(K) = (AN-0.1E+01)*(AN-0.2E+01)*PAR2
|
|
|
D1(K+1) = (AN+0.3E+01)*(AN+0.4E+01)*PAR2
|
|
|
V(K+3) = AS-(AN2-0.4E+01)*AC
|
|
|
AN = AN+0.2E+01
|
|
|
20 CONTINUE
|
|
|
AN2 = AN*AN
|
|
|
D(NOEQU) = -0.2E+01*(AN2-0.4E+01)*(PAR22-AN2-AN2)
|
|
|
V(NOEQU+3) = AS-(AN2-0.4E+01)*AC
|
|
|
V(4) = V(4)-0.56E+02*PAR2*V(3)
|
|
|
ASS = PARINT*SINPAR
|
|
|
ASAP = (((((0.210E+03*PAR2-0.1E+01)*COSPAR-(0.105E+03*PAR2
|
|
|
* -0.63E+02)*ASS)/AN2-(0.1E+01-0.15E+02*PAR2)*COSPAR
|
|
|
* +0.15E+02*ASS)/AN2-COSPAR+0.3E+01*ASS)/AN2-COSPAR)/AN2
|
|
|
V(NOEQU+3) = V(NOEQU+3)-0.2E+01*ASAP*PAR2*(AN-0.1E+01)*
|
|
|
* (AN-0.2E+01)
|
|
|
C
|
|
|
C SOLVE THE TRIDIAGONAL SYSTEM BY MEANS OF GAUSSIAN
|
|
|
C ELIMINATION WITH PARTIAL PIVOTING.
|
|
|
C
|
|
|
CALL SGTSL(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.4E+01
|
|
|
DO 40 I = 4,13
|
|
|
AN2 = AN*AN
|
|
|
V(I) = ((AN2-0.4E+01)*(0.2E+01*(PAR22-AN2-AN2)*V(I-1)-AC)
|
|
|
* +AS-PAR2*(AN+0.1E+01)*(AN+0.2E+01)*V(I-2))/
|
|
|
* (PAR2*(AN-0.1E+01)*(AN-0.2E+01))
|
|
|
AN = AN+0.2E+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.2E+01*(SINPAR-PARINT*COSPAR)/PAR2
|
|
|
V(2) = (0.18E+02-0.48E+02/PAR2)*SINPAR/PAR2
|
|
|
* +(-0.2E+01+0.48E+02/PAR2)*COSPAR/PARINT
|
|
|
AC = -0.24E+02*PARINT*COSPAR
|
|
|
AS = -0.8E+01*SINPAR
|
|
|
IF(ABS(PARINT).GT.0.24E+02) GO TO 80
|
|
|
C
|
|
|
C COMPUTE THE CHEBYSHEV MOMENTS AS THE
|
|
|
C SOLUTIONS OF A BOUNDARY VALUE PROBLEM WITH 1
|
|
|
C INITIAL VALUE (V(2)) AND 1 END VALUE (COMPUTED
|
|
|
C USING AN ASYMPTOTIC FORMULA).
|
|
|
C
|
|
|
AN = 0.5E+01
|
|
|
DO 70 K = 1,NOEQ1
|
|
|
AN2 = AN*AN
|
|
|
D(K) = -0.2E+01*(AN2-0.4E+01)*(PAR22-AN2-AN2)
|
|
|
D2(K) = (AN-0.1E+01)*(AN-0.2E+01)*PAR2
|
|
|
D1(K+1) = (AN+0.3E+01)*(AN+0.4E+01)*PAR2
|
|
|
V(K+2) = AC+(AN2-0.4E+01)*AS
|
|
|
AN = AN+0.2E+01
|
|
|
70 CONTINUE
|
|
|
AN2 = AN*AN
|
|
|
D(NOEQU) = -0.2E+01*(AN2-0.4E+01)*(PAR22-AN2-AN2)
|
|
|
V(NOEQU+2) = AC+(AN2-0.4E+01)*AS
|
|
|
V(3) = V(3)-0.42E+02*PAR2*V(2)
|
|
|
ASS = PARINT*COSPAR
|
|
|
ASAP = (((((0.105E+03*PAR2-0.63E+02)*ASS+(0.210E+03*PAR2
|
|
|
* -0.1E+01)*SINPAR)/AN2+(0.15E+02*PAR2-0.1E+01)*SINPAR-
|
|
|
* 0.15E+02*ASS)/AN2-0.3E+01*ASS-SINPAR)/AN2-SINPAR)/AN2
|
|
|
V(NOEQU+2) = V(NOEQU+2)-0.2E+01*ASAP*PAR2*(AN-0.1E+01)
|
|
|
* *(AN-0.2E+01)
|
|
|
C
|
|
|
C SOLVE THE TRIDIAGONAL SYSTEM BY MEANS OF GAUSSIAN
|
|
|
C ELIMINATION WITH PARTIAL PIVOTING.
|
|
|
C
|
|
|
CALL SGTSL(NOEQU,D1,D,D2,V(3),IERS)
|
|
|
GO TO 100
|
|
|
C
|
|
|
C COMPUTE THE CHEBYSHEV MOMENTS BY MEANS OF
|
|
|
C FORWARD RECURSION.
|
|
|
C
|
|
|
80 AN = 0.3E+01
|
|
|
DO 90 I = 3,12
|
|
|
AN2 = AN*AN
|
|
|
V(I) = ((AN2-0.4E+01)*(0.2E+01*(PAR22-AN2-AN2)*V(I-1)+AS)
|
|
|
* +AC-PAR2*(AN+0.1E+01)*(AN+0.2E+01)*V(I-2))
|
|
|
* /(PAR2*(AN-0.1E+01)*(AN-0.2E+01))
|
|
|
AN = AN+0.2E+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.5E+00*F(CENTR+HLGTH)
|
|
|
FVAL(13) = F(CENTR)
|
|
|
FVAL(25) = 0.5E+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 QCHEB(X,FVAL,CHEB12,CHEB24)
|
|
|
C
|
|
|
C COMPUTE THE INTEGRAL AND ERROR ESTIMATES.
|
|
|
C
|
|
|
RESC12 = CHEB12(13)*CHEBMO(M,13)
|
|
|
RESS12 = 0.0E+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.0E+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
|
|
|
c write(*,*) "1: ",RESULT
|
|
|
RESULT = CONC*RESC24-CONS*RESS24
|
|
|
ABSERR = ABS(CONC*ESTC)+ABS(CONS*ESTS)
|
|
|
GO TO 170
|
|
|
160 RESULT = CONC*RESS24+CONS*RESC24
|
|
|
c write(*,*) "2: ",RESULT
|
|
|
ABSERR = ABS(CONC*ESTS)+ABS(CONS*ESTC)
|
|
|
170 RETURN
|
|
|
END
|
|
|
|
