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