home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / fortran / library / ssp / numdiff / dbar.for next >
Text File  |  1986-01-03  |  4KB  |  131 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C     SUBROUTINE DBAR
  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 INTERVAL -
  10. C     THE SET OF T BETWEEN X AND X+H (H POSITIVE OR NEGATIVE) - USING
  11. C     FUNCTION VALUES ONLY ON THAT INTERVAL.
  12. C
  13. C      USAGE
  14. C     CALL DBAR(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 THAT DEFINES THE CLOSED INTERVAL WHOSE END-
  20. C           POINTS ARE X AND X+H (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 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 (MAGNITUDE OF THE) INTERNAL VALUE HH, WHICH IS DETER-
  32. C          MINED ACCORDING TO IH, IS THE MAXIMUM STEP-SIZE USED IN
  33. C          THE COMPUTATION OF THE ONE-SIDED DIVIDED DIFFERENCES (SEE
  34. C          METHOD.)    IF IH IS NON-ZERO, THEN THE SUBROUTINE GENERATES
  35. C          HH ACCORDING TO CRITERIA THAT BALANCE ROUND-OFF AND TRUN-
  36. C          CATION ERROR.  HH ALWAYS HAS THE SAME SIGN AS H AND IT IS
  37. C          ALWAYS LESS THAN OR EQUAL TO THE MAGNITUDE OF H IN AB-
  38. C          SOLUTE VALUE, SO THAT ALL COMPUTATION OCCURS IN THE CLOSED
  39. C          INTERVAL DETERMINED BY H.
  40. C
  41. C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  42. C     THE EXTERNAL FUNCTION SUBPROGRAM FCT(T) MUST BE FURNISHED BY
  43. C     THE USER.
  44. C
  45. C     METHOD
  46. C     THE COMPUTATION OF Z IS BASED ON RICHARDSON'S AND ROMBERG'S
  47. C     EXTRAPOLATION METHOD AS APPLIED TO THE SEQUENCE OF ONE-SIDED
  48. C     DIVIDED DIFFERENCES ASSOCIATED WITH THE POINT PAIRS
  49. C     (X,X+(K*HH)/10)K=1,...,10.  (SEE FILLIPI, S. AND ENGELS, H.,
  50. C     ALTES UND NEUES ZUR NUMERISCHEN DIFFERENTIATION, ELECTRONISCHE
  51. C     DATENVERARBEITUNG, ISS. 2 (1966), PP. 57-65.)
  52. C
  53. C     ..................................................................
  54. C
  55.       SUBROUTINE DBAR(X,H,IH,FCT,Z)
  56. C
  57. C
  58.       DIMENSION AUX(10)
  59. C
  60. C     NO ACTION IN CASE OF ZERO INTERVAL LENGTH
  61.       IF(H)1,17,1
  62. C
  63. C     GENERATE STEPSIZE HH FOR DIVIDED DIFFERENCES
  64.     1 C=ABS(H)
  65.       B=H
  66.       D=X
  67.       D=FCT(D)
  68.       IF(IH)2,9,2
  69.     2 HH=.5
  70.       IF(C-HH)3,4,4
  71.     3 HH=B
  72.     4 HH=SIGN(HH,B)
  73.       Z=ABS((FCT(X+HH)-D)/HH)
  74.       A=ABS(D)
  75.       HH=1.
  76.       IF(A-1.)6,6,5
  77.     5 HH=HH*A
  78.     6 IF(Z-1.)8,8,7
  79.     7 HH=HH/Z
  80.     8 IF(HH-C)10,10,9
  81.     9 HH=B
  82.    10 HH=SIGN(HH,B)
  83. C
  84. C     INITIALIZE DIFFERENTIATION LOOP
  85.       Z=(FCT(X+HH)-D)/HH
  86.       J=10
  87.       JJ=J-1
  88.       AUX(J)=Z
  89.       DH=HH/FLOAT(J)
  90.       DZ=1.E38
  91. C
  92. C     START DIFFERENTIATION LOOP
  93.    11 J=J-1
  94.       C=J
  95.       HH=C*DH
  96.       AUX(J)=(FCT(X+HH)-D)/HH
  97. C
  98. C     INITIALIZE EXTRAPOLATION LOOP
  99.       D2=1.E38
  100.       B=0.
  101.       A=1./C
  102. C
  103. C        START EXTRAPOLATION LOOP
  104.       DO 12 I=J,JJ
  105.       D1=D2
  106.       B=B+A
  107.       HH=(AUX(I)-AUX(I+1))/B
  108.       AUX(I+1)=AUX(I)+HH
  109. C
  110. C        TEST ON OSCILLATING INCREMENTS
  111.       D2=ABS(HH)
  112.       IF(D2-D1)12,13,13
  113.    12 CONTINUE
  114. C        END OF EXTRAPOLATION LOOP
  115. C
  116. C        UPDATE RESULT VALUE Z
  117.       I=JJ+1
  118.       GO TO 14
  119.    13 D2=D1
  120.       JJ=I
  121.    14 IF(D2-DZ)15,16,16
  122.    15 DZ=D2
  123.       Z=AUX(I)
  124.    16 IF(J-1)17,17,11
  125. C        END OF DIFFERENTIATION LOOP
  126. C
  127.    17 RETURN
  128.       END
  129. F EXTRAPOLATION LOOP
  130. C
  131. C        UPDATE RESULT VALUE