home *** CD-ROM | disk | FTP | other *** search
/ Frostbyte's 1980s DOS Shareware Collection / floppyshareware.zip / floppyshareware / DOOG / PCSSP1.ZIP / NUMQUAD.ZIP / DQATR.FOR next >
Text File  |  1985-11-29  |  4KB  |  120 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE DQATR
  5. C
  6. C
  7. C        PURPOSE
  8. C           TO COMPUTE AN APPROXIMATION FOR INTEGRAL(FCT(X), SUMMED
  9. C           OVER X FROM XL TO XU).
  10. C
  11. C        USAGE
  12. C           CALL DQATR (XL,XU,EPS,NDIM,FCT,Y,IER,AUX)
  13. C           PARAMETER FCT REQUIRES AN EXTERNAL STATEMENT.
  14. C
  15. C        DESCRIPTION OF PARAMETERS
  16. C           XL     - DOUBLE PRECISION LOWER BOUND OF THE INTERVAL.
  17. C           XU     - DOUBLE PRECISION UPPER BOUND OF THE INTERVAL.
  18. C           EPS    - SINGLE PRECISION UPPER BOUND OF THE ABSOLUTE ERROR.
  19. C           NDIM   - THE DIMENSION OF THE AUXILIARY STORAGE ARRAY AUX.
  20. C                    NDIM-1 IS THE MAXIMAL NUMBER OF BISECTIONS OF
  21. C                    THE INTERVAL (XL,XU).
  22. C           FCT    - THE NAME OF THE EXTERNAL DOUBLE PRECISION FUNCTION
  23. C                    SUBPROGRAM USED.
  24. C           Y      - RESULTING DOUBLE PRECISION APPROXIMATION FOR THE
  25. C                    INTEGRAL VALUE.
  26. C           IER    - A RESULTING ERROR PARAMETER.
  27. C           AUX    - AUXILIARY DOUBLE PRECISION STORAGE ARRAY WITH
  28. C                    DIMENSION NDIM.
  29. C
  30. C        REMARKS
  31. C           ERROR PARAMETER IER IS CODED IN THE FOLLOWING FORM
  32. C           IER=0  - IT WAS POSSIBLE TO REACH THE REQUIRED ACCURACY.
  33. C                    NO ERROR.
  34. C           IER=1  - IT IS IMPOSSIBLE TO REACH THE REQUIRED ACCURACY
  35. C                    BECAUSE OF ROUNDING ERRORS.
  36. C           IER=2  - IT WAS IMPOSSIBLE TO CHECK ACCURACY BECAUSE NDIM
  37. C                    IS LESS THAN 5, OR THE REQUIRED ACCURACY COULD NOT
  38. C                    BE REACHED WITHIN NDIM-1 STEPS. NDIM SHOULD BE
  39. C                    INCREASED.
  40. C
  41. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  42. C           THE EXTERNAL DOUBLE PRECISION FUNCTION SUBPROGRAM FCT(X)
  43. C           MUST BE CODED BY THE USER. ITS DOUBLE PRECISION ARGUMENT X
  44. C           SHOULD NOT BE DESTROYED.
  45. C
  46. C        METHOD
  47. C           EVALUATION OF Y IS DONE BY MEANS OF TRAPEZOIDAL RULE IN
  48. C           CONNECTION WITH ROMBERGS PRINCIPLE. ON RETURN Y CONTAINS
  49. C           THE BEST POSSIBLE APPROXIMATION OF THE INTEGRAL VALUE AND
  50. C           VECTOR AUX THE UPWARD DIAGONAL OF ROMBERG SCHEME.
  51. C           COMPONENTS AUX(I) (I=1,2,...,IEND, WITH IEND LESS THAN OR
  52. C           EQUAL TO NDIM) BECOME APPROXIMATIONS TO INTEGRAL VALUE WITH
  53. C           DECREASING ACCURACY BY MULTIPLICATION WITH (XU-XL).
  54. C           FOR REFERENCE, SEE
  55. C           (1) FILIPPI, DAS VERFAHREN VON ROMBERG-STIEFEL-BAUER ALS
  56. C               SPEZIALFALL DES ALLGEMEINEN PRINZIPS VON RICHARDSON,
  57. C               MATHEMATIK-TECHNIK-WIRTSCHAFT, VOL.11, ISS.2 (1964),
  58. C               PP.49-54.
  59. C           (2) BAUER, ALGORITHM 60, CACM, VOL.4, ISS.6 (1961), PP.255.
  60. C
  61. C     ..................................................................
  62. C
  63.       SUBROUTINE DQATR(XL,XU,EPS,NDIM,FCT,Y,IER,AUX)
  64. C
  65. C
  66.       DIMENSION AUX(1)
  67.       DOUBLE PRECISION AUX,XL,XU,X,Y,H,HH,HD,P,Q,SM,FCT
  68. C
  69. C     PREPARATIONS OF ROMBERG-LOOP
  70.       AUX(1)=.5D0*(FCT(XL)+FCT(XU))
  71.       H=XU-XL
  72.       IF(NDIM-1)8,8,1
  73.     1 IF(H)2,10,2
  74. C
  75. C     NDIM IS GREATER THAN 1 AND H IS NOT EQUAL TO 0.
  76.     2 HH=H
  77.       E=EPS/DABS(H)
  78.       DELT2=0.
  79.       P=1.D0
  80.       JJ=1
  81.       DO 7 I=2,NDIM
  82.       Y=AUX(1)
  83.       DELT1=DELT2
  84.       HD=HH
  85.       HH=.5D0*HH
  86.       P=.5D0*P
  87.       X=XL+HH
  88.       SM=0.D0
  89.       DO 3 J=1,JJ
  90.       SM=SM+FCT(X)
  91.     3 X=X+HD
  92.       AUX(I)=.5D0*AUX(I-1)+P*SM
  93. C     A NEW APPROXIMATION OF INTEGRAL VALUE IS COMPUTED BY MEANS OF
  94. C     TRAPEZOIDAL RULE.
  95. C
  96. C     START OF ROMBERGS EXTRAPOLATION METHOD.
  97.       Q=1.D0
  98.       JI=I-1
  99.       DO 4 J=1,JI
  100.       II=I-J
  101.       Q=Q+Q
  102.       Q=Q+Q
  103.     4 AUX(II)=AUX(II+1)+(AUX(II+1)-AUX(II))/(Q-1.D0)
  104. C     END OF ROMBERG-STEP
  105. C
  106.       DELT2=DABS(Y-AUX(1))
  107.       IF(I-5)7,5,5
  108.     5 IF(DELT2-E)10,10,6
  109.     6 IF(DELT2-DELT1)7,11,11
  110.     7 JJ=JJ+JJ
  111.     8 IER=2
  112.     9 Y=H*AUX(1)
  113.       RETURN
  114.    10 IER=0
  115.       GO TO 9
  116.    11 IER=1
  117.       Y=H*Y
  118.       RETURN
  119.       END
  120.