home *** CD-ROM | disk | FTP | other *** search
/ Frostbyte's 1980s DOS Shareware Collection / floppyshareware.zip / floppyshareware / DOOG / PCSSP1.ZIP / NUMQUAD.ZIP / DQSF.FOR < prev    next >
Text File  |  1985-11-29  |  3KB  |  118 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE DQSF
  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 DQSF (H,Y,Z,NDIM)
  12. C
  13. C        DESCRIPTION OF PARAMETERS
  14. C           H      - DOUBLE PRECISION INCREMENT OF ARGUMENT VALUES.
  15. C           Y      - DOUBLE PRECISION INPUT VECTOR OF FUNCTION VALUES.
  16. C           Z      - RESULTING DOUBLE PRECISION VECTOR OF INTEGRAL
  17. C                    VALUES. Z MAY BE 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 DQSF(H,Y,Z,NDIM)
  42. C
  43. C
  44.       DIMENSION Y(1),Z(1)
  45.       DOUBLE PRECISION Y,Z,H,HT,SUM1,SUM2,AUX,AUX1,AUX2
  46. C
  47.       HT=.33333333333333333D0*H
  48.       IF(NDIM-5)7,8,1
  49. C
  50. C     NDIM IS GREATER THAN 5. PREPARATIONS OF INTEGRATION LOOP
  51.     1 SUM1=Y(2)+Y(2)
  52.       SUM1=SUM1+SUM1
  53.       SUM1=HT*(Y(1)+SUM1+Y(3))
  54.       AUX1=Y(4)+Y(4)
  55.       AUX1=AUX1+AUX1
  56.       AUX1=SUM1+HT*(Y(3)+AUX1+Y(5))
  57.       AUX2=HT*(Y(1)+3.875D0*(Y(2)+Y(5))+2.625D0*(Y(3)+Y(4))+Y(6))
  58.       SUM2=Y(5)+Y(5)
  59.       SUM2=SUM2+SUM2
  60.       SUM2=AUX2-HT*(Y(4)+SUM2+Y(6))
  61.       Z(1)=0.D0
  62.       AUX=Y(3)+Y(3)
  63.       AUX=AUX+AUX
  64.       Z(2)=SUM2-HT*(Y(2)+AUX+Y(4))
  65.       Z(3)=SUM1
  66.       Z(4)=SUM2
  67.       IF(NDIM-6)5,5,2
  68. C
  69. C     INTEGRATION LOOP
  70.     2 DO 4 I=7,NDIM,2
  71.       SUM1=AUX1
  72.       SUM2=AUX2
  73.       AUX1=Y(I-1)+Y(I-1)
  74.       AUX1=AUX1+AUX1
  75.       AUX1=SUM1+HT*(Y(I-2)+AUX1+Y(I))
  76.       Z(I-2)=SUM1
  77.       IF(I-NDIM)3,6,6
  78.     3 AUX2=Y(I)+Y(I)
  79.       AUX2=AUX2+AUX2
  80.       AUX2=SUM2+HT*(Y(I-1)+AUX2+Y(I+1))
  81.     4 Z(I-1)=SUM2
  82.     5 Z(NDIM-1)=AUX1
  83.       Z(NDIM)=AUX2
  84.       RETURN
  85.     6 Z(NDIM-1)=SUM2
  86.       Z(NDIM)=AUX1
  87.       RETURN
  88. C     END OF INTEGRATION LOOP
  89. C
  90.     7 IF(NDIM-3)12,11,8
  91. C
  92. C     NDIM IS EQUAL TO 4 OR 5
  93.     8 SUM2=1.125D0*HT*(Y(1)+Y(2)+Y(2)+Y(2)+Y(3)+Y(3)+Y(3)+Y(4))
  94.       SUM1=Y(2)+Y(2)
  95.       SUM1=SUM1+SUM1
  96.       SUM1=HT*(Y(1)+SUM1+Y(3))
  97.       Z(1)=0.D0
  98.       AUX1=Y(3)+Y(3)
  99.       AUX1=AUX1+AUX1
  100.       Z(2)=SUM2-HT*(Y(2)+AUX1+Y(4))
  101.       IF(NDIM-5)10,9,9
  102.     9 AUX1=Y(4)+Y(4)
  103.       AUX1=AUX1+AUX1
  104.       Z(5)=SUM1+HT*(Y(3)+AUX1+Y(5))
  105.    10 Z(3)=SUM1
  106.       Z(4)=SUM2
  107.       RETURN
  108. C
  109. C     NDIM IS EQUAL TO 3
  110.    11 SUM1=HT*(1.25D0*Y(1)+Y(2)+Y(2)-.25D0*Y(3))
  111.       SUM2=Y(2)+Y(2)
  112.       SUM2=SUM2+SUM2
  113.       Z(3)=HT*(Y(1)+SUM2+Y(3))
  114.       Z(1)=0.D0
  115.       Z(2)=SUM1
  116.    12 RETURN
  117.       END
  118.