home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / quadpack / dqk21.f < prev    next >
Text File  |  1996-09-28  |  7KB  |  189 lines
  1.       SUBROUTINE DQK21(F,A,B,RESULT,ABSERR,RESABS,RESASC,IERR)
  2. C***BEGIN PROLOGUE  DQK21
  3. C***DATE WRITTEN   800101   (YYMMDD)
  4. C***REVISION DATE  830518   (YYMMDD)
  5. C***CATEGORY NO.  H2A1A2
  6. C***KEYWORDS  21-POINT GAUSS-KRONROD RULES
  7. C***AUTHOR  PIESSENS,ROBERT,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN
  8. C           DE DONCKER,ELISE,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN
  9. C***PURPOSE  TO COMPUTE I = INTEGRAL OF F OVER (A,B), WITH ERROR
  10. C                           ESTIMATE
  11. C                       J = INTEGRAL OF ABS(F) OVER (A,B)
  12. C***DESCRIPTION
  13. C
  14. C           INTEGRATION RULES
  15. C           STANDARD FORTRAN SUBROUTINE
  16. C           DOUBLE PRECISION VERSION
  17. C
  18. C           PARAMETERS
  19. C            ON ENTRY
  20. C              F      - DOUBLE PRECISION
  21. C                       FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
  22. C                       FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
  23. C                       DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
  24. C
  25. C              A      - DOUBLE PRECISION
  26. C                       LOWER LIMIT OF INTEGRATION
  27. C
  28. C              B      - DOUBLE PRECISION
  29. C                       UPPER LIMIT OF INTEGRATION
  30. C
  31. C            ON RETURN
  32. C              RESULT - DOUBLE PRECISION
  33. C                       APPROXIMATION TO THE INTEGRAL I
  34. C                       RESULT IS COMPUTED BY APPLYING THE 21-POINT
  35. C                       KRONROD RULE (RESK) OBTAINED BY OPTIMAL ADDITION
  36. C                       OF ABSCISSAE TO THE 10-POINT GAUSS RULE (RESG).
  37. C
  38. C              ABSERR - DOUBLE PRECISION
  39. C                       ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
  40. C                       WHICH SHOULD NOT EXCEED ABS(I-RESULT)
  41. C
  42. C              RESABS - DOUBLE PRECISION
  43. C                       APPROXIMATION TO THE INTEGRAL J
  44. C
  45. C              RESASC - DOUBLE PRECISION
  46. C                       APPROXIMATION TO THE INTEGRAL OF ABS(F-I/(B-A))
  47. C                       OVER (A,B)
  48. C
  49. C***REFERENCES  (NONE)
  50. C***ROUTINES CALLED  D1MACH
  51. C***END PROLOGUE  DQK21
  52. C
  53.       DOUBLE PRECISION A,ABSC,ABSERR,B,CENTR,DABS,DHLGTH,DMAX1,DMIN1,
  54.      *  D1MACH,EPMACH,F,FC,FSUM,FVAL1,FVAL2,FV1,FV2,HLGTH,RESABS,RESASC,
  55.      *  RESG,RESK,RESKH,RESULT,UFLOW,WG,WGK,XGK
  56.       INTEGER J,JTW,JTWM1
  57.       EXTERNAL F
  58. C
  59.       DIMENSION FV1(10),FV2(10),WG(5),WGK(11),XGK(11)
  60. C
  61. C           THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE INTERVAL (-1,1).
  62. C           BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND THEIR
  63. C           CORRESPONDING WEIGHTS ARE GIVEN.
  64. C
  65. C           XGK    - ABSCISSAE OF THE 21-POINT KRONROD RULE
  66. C                    XGK(2), XGK(4), ...  ABSCISSAE OF THE 10-POINT
  67. C                    GAUSS RULE
  68. C                    XGK(1), XGK(3), ...  ABSCISSAE WHICH ARE OPTIMALLY
  69. C                    ADDED TO THE 10-POINT GAUSS RULE
  70. C
  71. C           WGK    - WEIGHTS OF THE 21-POINT KRONROD RULE
  72. C
  73. C           WG     - WEIGHTS OF THE 10-POINT GAUSS RULE
  74. C
  75. C
  76. C GAUSS QUADRATURE WEIGHTS AND KRONRON QUADRATURE ABSCISSAE AND WEIGHTS
  77. C AS EVALUATED WITH 80 DECIMAL DIGIT ARITHMETIC BY L. W. FULLERTON,
  78. C BELL LABS, NOV. 1981.
  79. C
  80.       DATA WG  (  1) / 0.0666713443 0868813759 3568809893 332 D0 /
  81.       DATA WG  (  2) / 0.1494513491 5058059314 5776339657 697 D0 /
  82.       DATA WG  (  3) / 0.2190863625 1598204399 5534934228 163 D0 /
  83.       DATA WG  (  4) / 0.2692667193 0999635509 1226921569 469 D0 /
  84.       DATA WG  (  5) / 0.2955242247 1475287017 3892994651 338 D0 /
  85. C
  86.       DATA XGK (  1) / 0.9956571630 2580808073 5527280689 003 D0 /
  87.       DATA XGK (  2) / 0.9739065285 1717172007 7964012084 452 D0 /
  88.       DATA XGK (  3) / 0.9301574913 5570822600 1207180059 508 D0 /
  89.       DATA XGK (  4) / 0.8650633666 8898451073 2096688423 493 D0 /
  90.       DATA XGK (  5) / 0.7808177265 8641689706 3717578345 042 D0 /
  91.       DATA XGK (  6) / 0.6794095682 9902440623 4327365114 874 D0 /
  92.       DATA XGK (  7) / 0.5627571346 6860468333 9000099272 694 D0 /
  93.       DATA XGK (  8) / 0.4333953941 2924719079 9265943165 784 D0 /
  94.       DATA XGK (  9) / 0.2943928627 0146019813 1126603103 866 D0 /
  95.       DATA XGK ( 10) / 0.1488743389 8163121088 4826001129 720 D0 /
  96.       DATA XGK ( 11) / 0.0000000000 0000000000 0000000000 000 D0 /
  97. C
  98.       DATA WGK (  1) / 0.0116946388 6737187427 8064396062 192 D0 /
  99.       DATA WGK (  2) / 0.0325581623 0796472747 8818972459 390 D0 /
  100.       DATA WGK (  3) / 0.0547558965 7435199603 1381300244 580 D0 /
  101.       DATA WGK (  4) / 0.0750396748 1091995276 7043140916 190 D0 /
  102.       DATA WGK (  5) / 0.0931254545 8369760553 5065465083 366 D0 /
  103.       DATA WGK (  6) / 0.1093871588 0229764189 9210590325 805 D0 /
  104.       DATA WGK (  7) / 0.1234919762 6206585107 7958109831 074 D0 /
  105.       DATA WGK (  8) / 0.1347092173 1147332592 8054001771 707 D0 /
  106.       DATA WGK (  9) / 0.1427759385 7706008079 7094273138 717 D0 /
  107.       DATA WGK ( 10) / 0.1477391049 0133849137 4841515972 068 D0 /
  108.       DATA WGK ( 11) / 0.1494455540 0291690566 4936468389 821 D0 /
  109. C
  110. C
  111. C           LIST OF MAJOR VARIABLES
  112. C           -----------------------
  113. C
  114. C           CENTR  - MID POINT OF THE INTERVAL
  115. C           HLGTH  - HALF-LENGTH OF THE INTERVAL
  116. C           ABSC   - ABSCISSA
  117. C           FVAL*  - FUNCTION VALUE
  118. C           RESG   - RESULT OF THE 10-POINT GAUSS FORMULA
  119. C           RESK   - RESULT OF THE 21-POINT KRONROD FORMULA
  120. C           RESKH  - APPROXIMATION TO THE MEAN VALUE OF F OVER (A,B),
  121. C                    I.E. TO I/(B-A)
  122. C
  123. C
  124. C           MACHINE DEPENDENT CONSTANTS
  125. C           ---------------------------
  126. C
  127. C           EPMACH IS THE LARGEST RELATIVE SPACING.
  128. C           UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
  129. C
  130. C***FIRST EXECUTABLE STATEMENT  DQK21
  131.       EPMACH = D1MACH(4)
  132.       UFLOW = D1MACH(1)
  133. C
  134.       CENTR = 0.5D+00*(A+B)
  135.       HLGTH = 0.5D+00*(B-A)
  136.       DHLGTH = DABS(HLGTH)
  137. C
  138. C           COMPUTE THE 21-POINT KRONROD APPROXIMATION TO
  139. C           THE INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR.
  140. C
  141.       RESG = 0.0D+00
  142.       IERR = 0
  143.       FC = F(CENTR,IERR)
  144.       IF (IERR .LT. 0) RETURN
  145.       RESK = WGK(11)*FC
  146.       RESABS = DABS(RESK)
  147.       DO 10 J=1,5
  148.         JTW = 2*J
  149.         ABSC = HLGTH*XGK(JTW)
  150.         FVAL1 = F(CENTR-ABSC,IERR)
  151.         IF (IERR .LT. 0) RETURN
  152.         FVAL2 = F(CENTR+ABSC,IERR)
  153.         IF (IERR .LT. 0) RETURN
  154.         FV1(JTW) = FVAL1
  155.         FV2(JTW) = FVAL2
  156.         FSUM = FVAL1+FVAL2
  157.         RESG = RESG+WG(J)*FSUM
  158.         RESK = RESK+WGK(JTW)*FSUM
  159.         RESABS = RESABS+WGK(JTW)*(DABS(FVAL1)+DABS(FVAL2))
  160.    10 CONTINUE
  161.       DO 15 J = 1,5
  162.         JTWM1 = 2*J-1
  163.         ABSC = HLGTH*XGK(JTWM1)
  164.         FVAL1 = F(CENTR-ABSC,IERR)
  165.         IF (IERR .LT. 0) RETURN
  166.         FVAL2 = F(CENTR+ABSC,IERR)
  167.         IF (IERR .LT. 0) RETURN
  168.         FV1(JTWM1) = FVAL1
  169.         FV2(JTWM1) = FVAL2
  170.         FSUM = FVAL1+FVAL2
  171.         RESK = RESK+WGK(JTWM1)*FSUM
  172.         RESABS = RESABS+WGK(JTWM1)*(DABS(FVAL1)+DABS(FVAL2))
  173.    15 CONTINUE
  174.       RESKH = RESK*0.5D+00
  175.       RESASC = WGK(11)*DABS(FC-RESKH)
  176.       DO 20 J=1,10
  177.         RESASC = RESASC+WGK(J)*(DABS(FV1(J)-RESKH)+DABS(FV2(J)-RESKH))
  178.    20 CONTINUE
  179.       RESULT = RESK*HLGTH
  180.       RESABS = RESABS*DHLGTH
  181.       RESASC = RESASC*DHLGTH
  182.       ABSERR = DABS((RESK-RESG)*HLGTH)
  183.       IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.0D+00)
  184.      *  ABSERR = RESASC*DMIN1(0.1D+01,(0.2D+03*ABSERR/RESASC)**1.5D+00)
  185.       IF(RESABS.GT.UFLOW/(0.5D+02*EPMACH)) ABSERR = DMAX1
  186.      *  ((EPMACH*0.5D+02)*RESABS,ABSERR)
  187.       RETURN
  188.       END
  189.