Logo Search packages:      
Sourcecode: maxima version File versions  Download package

dqc25c.f

*DECK DQC25C
      SUBROUTINE DQC25C (F, A, B, C, RESULT, ABSERR, KRUL, NEVAL)
C***BEGIN PROLOGUE  DQC25C
C***PURPOSE  To compute I = Integral of F*W over (A,B) with
C            error estimate, where W(X) = 1/(X-C)
C***LIBRARY   SLATEC (QUADPACK)
C***CATEGORY  H2A2A2, J4
C***TYPE      DOUBLE PRECISION (QC25C-S, DQC25C-D)
C***KEYWORDS  25-POINT CLENSHAW-CURTIS INTEGRATION, 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 the computation of CAUCHY
C        PRINCIPAL VALUE integrals
C        Standard fortran subroutine
C        Double precision version
C
C        PARAMETERS
C           F      - Double precision
C                    Function subprogram defining the integrand function
C                    F(X). The actual name for F needs to be declared
C                    E X T E R N A L  in the driver program.
C
C           A      - Double precision
C                    Left end point of the integration interval
C
C           B      - Double precision
C                    Right end point of the integration interval, B.GT.A
C
C           C      - Double precision
C                    Parameter in the WEIGHT function
C
C           RESULT - Double precision
C                    Approximation to the integral
C                    result is computed by using a generalized
C                    Clenshaw-Curtis method if C lies within ten percent
C                    of the integration interval. In the other case the
C                    15-point Kronrod rule obtained by optimal addition
C                    of abscissae to the 7-point Gauss rule, is applied.
C
C           ABSERR - Double precision
C                    Estimate of the modulus of the absolute error,
C                    which should equal or exceed ABS(I-RESULT)
C
C           KRUL   - Integer
C                    Key which is decreased by 1 if the 15-point
C                    Gauss-Kronrod scheme has been used
C
C           NEVAL  - Integer
C                    Number of integrand evaluations
C
C ......................................................................
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  DQCHEB, DQK15W, DQWGTC
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  DQC25C
C
      DOUBLE PRECISION A,ABSERR,AK22,AMOM0,AMOM1,AMOM2,B,C,CC,CENTR,
     1  CHEB12,CHEB24,DQWGTC,F,FVAL,HLGTH,P2,P3,P4,RESABS,
     2  RESASC,RESULT,RES12,RES24,U,X
      INTEGER I,ISYM,K,KP,KRUL,NEVAL
C
      DIMENSION X(11),FVAL(25),CHEB12(13),CHEB24(25)
C
      EXTERNAL F, DQWGTC
C
C           THE VECTOR X CONTAINS THE VALUES COS(K*PI/24),
C           K = 1, ..., 11, TO BE USED FOR THE CHEBYSHEV SERIES
C           EXPANSION OF F
C
      SAVE X
      DATA X(1) / 0.9914448613 7381041114 4557526928 563D0 /
      DATA X(2) / 0.9659258262 8906828674 9743199728 897D0 /
      DATA X(3) / 0.9238795325 1128675612 8183189396 788D0 /
      DATA X(4) / 0.8660254037 8443864676 3723170752 936D0 /
      DATA X(5) / 0.7933533402 9123516457 9776961501 299D0 /
      DATA X(6) / 0.7071067811 8654752440 0844362104 849D0 /
      DATA X(7) / 0.6087614290 0872063941 6097542898 164D0 /
      DATA X(8) / 0.5000000000 0000000000 0000000000 000D0 /
      DATA X(9) / 0.3826834323 6508977172 8459984030 399D0 /
      DATA X(10) / 0.2588190451 0252076234 8898837624 048D0 /
      DATA X(11) / 0.1305261922 2005159154 8406227895 489D0 /
C
C           LIST OF MAJOR VARIABLES
C           ----------------------
C           FVAL   - VALUE OF THE FUNCTION F AT THE POINTS
C                    COS(K*PI/24),  K = 0, ..., 24
C           CHEB12 - CHEBYSHEV SERIES EXPANSION COEFFICIENTS,
C                    FOR THE FUNCTION F, OF DEGREE 12
C           CHEB24 - CHEBYSHEV SERIES EXPANSION COEFFICIENTS,
C                    FOR THE FUNCTION F, OF DEGREE 24
C           RES12  - APPROXIMATION TO THE INTEGRAL CORRESPONDING
C                    TO THE USE OF CHEB12
C           RES24  - APPROXIMATION TO THE INTEGRAL CORRESPONDING
C                    TO THE USE OF CHEB24
C           DQWGTC - EXTERNAL FUNCTION SUBPROGRAM DEFINING
C                    THE WEIGHT FUNCTION
C           HLGTH  - HALF-LENGTH OF THE INTERVAL
C           CENTR  - MID POINT OF THE INTERVAL
C
C
C           CHECK THE POSITION OF C.
C
C***FIRST EXECUTABLE STATEMENT  DQC25C
      CC = (0.2D+01*C-B-A)/(B-A)
      IF(ABS(CC).LT.0.11D+01) GO TO 10
C
C           APPLY THE 15-POINT GAUSS-KRONROD SCHEME.
C
      KRUL = KRUL-1
      CALL DQK15W(F,DQWGTC,C,P2,P3,P4,KP,A,B,RESULT,ABSERR,
     1  RESABS,RESASC)
      NEVAL = 15
      IF (RESASC.EQ.ABSERR) KRUL = KRUL+1
      GO TO 50
C
C           USE THE GENERALIZED CLENSHAW-CURTIS METHOD.
C
   10 HLGTH = 0.5D+00*(B-A)
      CENTR = 0.5D+00*(B+A)
      NEVAL = 25
      FVAL(1) = 0.5D+00*F(HLGTH+CENTR)
      FVAL(13) = F(CENTR)
      FVAL(25) = 0.5D+00*F(CENTR-HLGTH)
      DO 20 I=2,12
        U = HLGTH*X(I-1)
        ISYM = 26-I
        FVAL(I) = F(U+CENTR)
        FVAL(ISYM) = F(CENTR-U)
   20 CONTINUE
C
C           COMPUTE THE CHEBYSHEV SERIES EXPANSION.
C
      CALL DQCHEB(X,FVAL,CHEB12,CHEB24)
C
C           THE MODIFIED CHEBYSHEV MOMENTS ARE COMPUTED BY FORWARD
C           RECURSION, USING AMOM0 AND AMOM1 AS STARTING VALUES.
C
      AMOM0 = LOG(ABS((0.1D+01-CC)/(0.1D+01+CC)))
      AMOM1 = 0.2D+01+CC*AMOM0
      RES12 = CHEB12(1)*AMOM0+CHEB12(2)*AMOM1
      RES24 = CHEB24(1)*AMOM0+CHEB24(2)*AMOM1
      DO 30 K=3,13
        AMOM2 = 0.2D+01*CC*AMOM1-AMOM0
        AK22 = (K-2)*(K-2)
        IF((K/2)*2.EQ.K) AMOM2 = AMOM2-0.4D+01/(AK22-0.1D+01)
        RES12 = RES12+CHEB12(K)*AMOM2
        RES24 = RES24+CHEB24(K)*AMOM2
        AMOM0 = AMOM1
        AMOM1 = AMOM2
   30 CONTINUE
      DO 40 K=14,25
        AMOM2 = 0.2D+01*CC*AMOM1-AMOM0
        AK22 = (K-2)*(K-2)
        IF((K/2)*2.EQ.K) AMOM2 = AMOM2-0.4D+01/(AK22-0.1D+01)
        RES24 = RES24+CHEB24(K)*AMOM2
        AMOM0 = AMOM1
        AMOM1 = AMOM2
   40 CONTINUE
      RESULT = RES24
      ABSERR = ABS(RES24-RES12)
   50 RETURN
      END

Generated by  Doxygen 1.6.0   Back to index