home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / fortran / library / ssp / numquad / qsf.for < prev    next >
Encoding:
Text File  |  1985-11-29  |  3.2 KB  |  117 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE QSF
  5. C
  6. C        PURPOSE
  7. C           TO COMPUTE THE VECTOR OF INTEGRAL VALUES FOR A GIVEN
  8. C           EQUIDISTANT TABLE OF FUNCTION VALUES.
  9. C
  10. C        USAGE
  11. C           CALL QSF (H,Y,Z,NDIM)
  12. C
  13. C        DESCRIPTION OF PARAMETERS
  14. C           H      - THE INCREMENT OF ARGUMENT VALUES.
  15. C           Y      - THE INPUT VECTOR OF FUNCTION VALUES.
  16. C           Z      - THE RESULTING VECTOR OF INTEGRAL VALUES. Z MAY BE
  17. C                    IDENTICAL WITH Y.
  18. C           NDIM   - THE DIMENSION OF VECTORS Y AND Z.
  19. C
  20. C        REMARKS
  21. C           NO ACTION IN CASE NDIM LESS THAN 3.
  22. C
  23. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  24. C           NONE
  25. C
  26. C        METHOD
  27. C           BEGINNING WITH Z(1)=0, EVALUATION OF VECTOR Z IS DONE BY
  28. C           MEANS OF SIMPSONS RULE TOGETHER WITH NEWTONS 3/8 RULE OR A
  29. C           COMBINATION OF THESE TWO RULES. TRUNCATION ERROR IS OF
  30. C           ORDER H**5 (I.E. FOURTH ORDER METHOD). ONLY IN CASE NDIM=3
  31. C           TRUNCATION ERROR OF Z(2) IS OF ORDER H**4.
  32. C           FOR REFERENCE, SEE
  33. C           (1) F.B.HILDEBRAND, INTRODUCTION TO NUMERICAL ANALYSIS,
  34. C               MCGRAW-HILL, NEW YORK/TORONTO/LONDON, 1956, PP.71-76.
  35. C           (2) R.ZURMUEHL, PRAKTISCHE MATHEMATIK FUER INGENIEURE UND
  36. C               PHYSIKER, SPRINGER, BERLIN/GOETTINGEN/HEIDELBERG, 1963,
  37. C               PP.214-221.
  38. C
  39. C     ..................................................................
  40. C
  41.       SUBROUTINE QSF(H,Y,Z,NDIM)
  42. C
  43. C
  44.       DIMENSION Y(1),Z(1)
  45. C
  46.       HT=.3333333*H
  47.       IF(NDIM-5)7,8,1
  48. C
  49. C     NDIM IS GREATER THAN 5. PREPARATIONS OF INTEGRATION LOOP
  50.     1 SUM1=Y(2)+Y(2)
  51.       SUM1=SUM1+SUM1
  52.       SUM1=HT*(Y(1)+SUM1+Y(3))
  53.       AUX1=Y(4)+Y(4)
  54.       AUX1=AUX1+AUX1
  55.       AUX1=SUM1+HT*(Y(3)+AUX1+Y(5))
  56.       AUX2=HT*(Y(1)+3.875*(Y(2)+Y(5))+2.625*(Y(3)+Y(4))+Y(6))
  57.       SUM2=Y(5)+Y(5)
  58.       SUM2=SUM2+SUM2
  59.       SUM2=AUX2-HT*(Y(4)+SUM2+Y(6))
  60.       Z(1)=0.
  61.       AUX=Y(3)+Y(3)
  62.       AUX=AUX+AUX
  63.       Z(2)=SUM2-HT*(Y(2)+AUX+Y(4))
  64.       Z(3)=SUM1
  65.       Z(4)=SUM2
  66.       IF(NDIM-6)5,5,2
  67. C
  68. C     INTEGRATION LOOP
  69.     2 DO 4 I=7,NDIM,2
  70.       SUM1=AUX1
  71.       SUM2=AUX2
  72.       AUX1=Y(I-1)+Y(I-1)
  73.       AUX1=AUX1+AUX1
  74.       AUX1=SUM1+HT*(Y(I-2)+AUX1+Y(I))
  75.       Z(I-2)=SUM1
  76.       IF(I-NDIM)3,6,6
  77.     3 AUX2=Y(I)+Y(I)
  78.       AUX2=AUX2+AUX2
  79.       AUX2=SUM2+HT*(Y(I-1)+AUX2+Y(I+1))
  80.     4 Z(I-1)=SUM2
  81.     5 Z(NDIM-1)=AUX1
  82.       Z(NDIM)=AUX2
  83.       RETURN
  84.     6 Z(NDIM-1)=SUM2
  85.       Z(NDIM)=AUX1
  86.       RETURN
  87. C     END OF INTEGRATION LOOP
  88. C
  89.     7 IF(NDIM-3)12,11,8
  90. C
  91. C     NDIM IS EQUAL TO 4 OR 5
  92.     8 SUM2=1.125*HT*(Y(1)+Y(2)+Y(2)+Y(2)+Y(3)+Y(3)+Y(3)+Y(4))
  93.       SUM1=Y(2)+Y(2)
  94.       SUM1=SUM1+SUM1
  95.       SUM1=HT*(Y(1)+SUM1+Y(3))
  96.       Z(1)=0.
  97.       AUX1=Y(3)+Y(3)
  98.       AUX1=AUX1+AUX1
  99.       Z(2)=SUM2-HT*(Y(2)+AUX1+Y(4))
  100.       IF(NDIM-5)10,9,9
  101.     9 AUX1=Y(4)+Y(4)
  102.       AUX1=AUX1+AUX1
  103.       Z(5)=SUM1+HT*(Y(3)+AUX1+Y(5))
  104.    10 Z(3)=SUM1
  105.       Z(4)=SUM2
  106.       RETURN
  107. C
  108. C     NDIM IS EQUAL TO 3
  109.    11 SUM1=HT*(1.25*Y(1)+Y(2)+Y(2)-.25*Y(3))
  110.       SUM2=Y(2)+Y(2)
  111.       SUM2=SUM2+SUM2
  112.       Z(3)=HT*(Y(1)+SUM2+Y(3))
  113.       Z(1)=0.
  114.       Z(2)=SUM1
  115.    12 RETURN
  116.       END
  117.