home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / fortran / library / ssp / numdiff / dcar.for < prev    next >
Text File  |  1986-01-03  |  4KB  |  132 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C     SUBROUTINE DCAR
  5. C
  6. C     PURPOSE
  7. C     TO COMPUTE, AT A GIVEN POINT X, AN APPROXIMATION Z TO THE
  8. C     DERIVATIVE OF AN ANALYTICALLY GIVEN FUNCTION FCT THAT IS 11-
  9. C     TIMES DIFFERENTIABLE IN A DOMAIN CONTAINING A CLOSED, 2-SIDED
  10. C     SYMMETRIC INTERVAL OF RADIUS ABSOLUTE H ABOUT X, USING FUNCTION
  11. C     VALUES ONLY ON THAT CLOSED INTERVAL.
  12. C
  13. C     USAGE
  14. C        CALL DCAR (X,H,IH,FCT,Z)
  15. C     PARAMETER FCT REQUIRES AN EXTERNAL STATEMENT
  16. C
  17. C     DESCRIPTION OF PARAMETERS
  18. C     X   - THE POINT AT WHICH THE DERIVATIVE IS TO BE COMPUTED
  19. C     H   - THE NUMBER WHOSE ABSOLUTE VALUE DEFINES THE CLOSED,
  20. C           SYMMETRIC 2-SIDED INTERVAL ABOUT X (SEE PURPOSE)
  21. C     IH  - INPUT PARAMETER (SEE REMARKS AND METHOD)
  22. C           IH NON-ZERO - THE SUBROUTINE GENERATES THE INTERNAL
  23. C                 VALUE HH
  24. C           IH    =     0 - THE INTERNAL VALUE HH IS SET TO ABSOLUTE H
  25. C     FCT - THE NAME OF THE EXTERNAL FUNCTION SUBPROGRAM THAT WILL
  26. C           GENERATE THE NECESSARY FUNCTION VALUES
  27. C     Z   - RESULTING DERIVATIVE VALUE
  28. C
  29. C     REMARKS
  30. C     (1)  IF H = 0, THEN THERE IS NO COMPUTATION.
  31. C     (2)  THE INTERNAL VALUE HH, WHICH IS DETERMINED ACCORDING TO
  32. C          IH, IS THE MAXIMUM STEP-SIZE USED IN THE COMPUTATION OF
  33. C          THE CENTRAL DIVIDED DIFFERENCES (SEE METHOD.)  IF IH IS
  34. C          NON-ZERO, THEN THE SUBROUTINE GENERATES HH ACCORDING TO
  35. C          CRITERIA THAT BALANCE ROUND-OFF AND TRUNCATION ERROR.  HH
  36. C          IS ALWAYS LESS THAN OR EQUAL TO ABSOLUTE H IN ABSOLUTE
  37. C          VALUE, SO THAT ALL COMPUTATION OCCURS WITHIN A RADIUS
  38. C          ABSOLUTE H OF X.
  39. C
  40. C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  41. C     THE EXTERNAL FUNCTION SUBPROGRAM FCT(T) MUST BE FURNISHED BY
  42. C     THE USER.
  43. C
  44. C     METHOD
  45. C     THE COMPUTATION OF Z IS BASED ON RICHARDSON'S AND ROMBERG'S
  46. C     EXTRAPOLATION METHOD AS APPLIED TO THE SEQUENCE OF CENTRAL
  47. C     DIVIDED DIFFERENCES ASSOCIATED WITH THE POINT PAIRS
  48. C     (X-(K*HH)/5,X+(K*HH)/5) K=1,...,5.  (SEE FILLIPI, S. AND
  49. C     ENGELS, H., ALTES UND NEUES ZUR NUMERISCHEN DIFFERENTIATION,
  50. C     ELECTRONISCHE DATENVERARBEITUNG, ISS. 2 (1966), PP. 57-65.)
  51. C
  52. C     ..................................................................
  53. C
  54.       SUBROUTINE DCAR(X,H,IH,FCT,Z)
  55. C
  56. C
  57.       DIMENSION AUX(5)
  58. C
  59. C     NO ACTION IN CASE OF ZERO INTERVAL LENGTH
  60.       IF(H)1,17,1
  61. C
  62. C     GENERATE STEPSIZE HH FOR DIVIDED DIFFERENCES
  63.     1 C=ABS(H)
  64.       IF(IH)2,9,2
  65.     2 HH=.5
  66.       IF(C-HH)3,4,4
  67.     3 HH=C
  68.     4 A=FCT(X+HH)
  69.       B=FCT(X-HH)
  70.       Z=ABS((A-B)/(HH+HH))
  71.       A=.5*ABS(A+B)
  72.       HH=.5
  73.       IF(A-1.)6,6,5
  74.     5 HH=HH*A
  75.     6 IF(Z-1.)8,8,7
  76.     7 HH=HH/Z
  77.     8 IF(HH-C)10,10,9
  78.     9 HH=C
  79. C
  80. C     INITIALIZE DIFFERENTIATION LOOP
  81.    10 Z=(FCT(X+HH)-FCT(X-HH))/(HH+HH)
  82.       J=5
  83.       JJ=J-1
  84.       AUX(J)=Z
  85.       DH=HH/FLOAT(J)
  86.       DZ=1.E38
  87. C
  88. C     START DIFFERENTIATION LOOP
  89.    11 J=J-1
  90.       C=J
  91.       HH=C*DH
  92.       AUX(J)=(FCT(X+HH)-FCT(X-HH))/(HH+HH)
  93. C
  94. C     INITIALIZE EXTRAPOLATION LOOP
  95.       D2=1.E38
  96.       B=0.
  97.       A=1./C
  98. C
  99. C        START EXTRAPOLATION LOOP
  100.       DO 12 I=J,JJ
  101.       D1=D2
  102.       B=B+A
  103.       HH=(AUX(I)-AUX(I+1))/(B*(2.+B))
  104.       AUX(I+1)=AUX(I)+HH
  105. C
  106. C        TEST ON OSCILLATING INCREMENTS
  107.       D2=ABS(HH)
  108.       IF(D2-D1)12,13,13
  109.    12 CONTINUE
  110. C        END OF EXTRAPOLATION LOOP
  111. C
  112. C        UPDATE RESULT VALUE Z
  113.       I=JJ+1
  114.       GO TO 14
  115.    13 D2=D1
  116.       JJ=I
  117.    14 IF(D2-DZ)15,16,16
  118.    15 DZ=D2
  119.       Z=AUX(I)
  120.    16 IF(J-1)17,17,11
  121. C        END OF DIFFERENTIATION LOOP
  122. C
  123.    17 RETURN
  124.       END
  125.  EXTRAPOLATION LOOP
  126. C
  127. C        UPDATE RESULT VALUE Z
  128.       I=JJ+1
  129.       GO TO 14
  130.    13 D2=D1
  131.       JJ=I
  132.    14 IF(D2-DZ