home *** CD-ROM | disk | FTP | other *** search
/ Frostbyte's 1980s DOS Shareware Collection / floppyshareware.zip / floppyshareware / DOOG / PCSSP1.ZIP / NUMDIFF.ZIP / DDCAR.FOR < prev    next >
Text File  |  1986-01-03  |  4KB  |  130 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C     SUBROUTINE DDCAR
  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 DDCAR(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           X IS IN DOUBLE PRECISION.
  20. C     H   - THE NUMBER WHOSE ABSOLUTE VALUE DEFINES THE CLOSED,
  21. C           SYMMETRIC 2-SIDED INTERVAL ABOUT X (SEE PURPOSE)
  22. C           H IS IN SINGLE PRECISION
  23. C     IH  - INPUT PARAMETER (SEE REMARKS AND METHOD)
  24. C           IH NON-ZERO - THE SUBROUTINE GENERATES THE INTERNAL
  25. C                 VALUE HH
  26. C           IH    =     0 - THE INTERNAL VALUE HH IS SET TO ABSOLUTE H
  27. C     FCT - THE NAME OF THE EXTERNAL DOUBLE PRECISION FUNCTION
  28. C           SUBPROGRAM THAT WILL GENERATE THE NECESSARY FUNCTION
  29. C           VALUES.
  30. C     Z   - RESULTING DERIVATIVE VALUE - DOUBLE PRECISION
  31. C
  32. C     REMARKS
  33. C     (1)  IF H = 0, THEN THERE IS NO COMPUTATION.
  34. C     (2)  THE INTERNAL VALUE HH, WHICH IS DETERMINED ACCORDING TO
  35. C          IH, IS THE MAXIMUM STEP-SIZE USED IN THE COMPUTATION OF
  36. C          THE CENTRAL DIVIDED DIFFERENCES (SEE METHOD.)  IF IH IS
  37. C          NON-ZERO, THEN THE SUBROUTINE GENERATES HH ACCORDING TO
  38. C          CRITERIA THAT BALANCE ROUND-OFF AND TRUNCATION ERROR.  HH
  39. C          IS ALWAYS LESS THAN OR EQUAL TO ABSOLUTE H IN ABSOLUTE
  40. C          VALUE, SO THAT ALL COMPUTATION OCCURS WITHIN A RADIUS
  41. C          ABSOLUTE H OF X.
  42. C
  43. C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  44. C     THE EXTERNAL FUNCTION SUBPROGRAM FCT(T) MUST BE FURNISHED BY
  45. C     THE USER. FCT(T) IS IN DOUBLE PRECISION
  46. C
  47. C     METHOD
  48. C     THE COMPUTATION OF Z IS BASED ON RICHARDSON'S AND ROMBERG'S
  49. C     EXTRAPOLATION METHOD AS APPLIED TO THE SEQUENCE OF CENTRAL
  50. C     DIVIDED DIFFERENCES ASSOCIATED WITH THE POINT PAIRS
  51. C     (X-(K*HH)/5,X+(K*HH)/5) K=1,...,5.  (SEE FILLIPI, S. AND
  52. C     ENGELS, H., ALTES UND NEUES ZUR NUMERISCHEN DIFFERENTIATION,
  53. C     ELECTRONISCHE DATENVERARBEITUNG, ISS. 2 (1966), PP. 57-65.)
  54. C
  55. C     ..................................................................
  56. C
  57.       SUBROUTINE DDCAR(X,H,IH,FCT,Z)
  58. C
  59. C
  60.       DIMENSION AUX(5)
  61.       DOUBLE PRECISION X,FCT,Z,AUX,A,B,C,DH,HH
  62. C
  63. C     NO ACTION IN CASE OF ZERO INTERVAL LENGTH
  64.       IF(H)1,17,1
  65. C
  66. C     GENERATE STEPSIZE HH FOR DIVIDED DIFFERENCES
  67.     1 C=ABS(H)
  68.       IF(IH)2,9,2
  69.     2 HH=.5D-2
  70.       IF(C-HH)3,4,4
  71.     3 HH=C
  72.     4 A=FCT(X+HH)
  73.       B=FCT(X-HH)
  74.       Z=DABS((A-B)/(HH+HH))
  75.       A=.5D0*DABS(A+B)
  76.       HH=.5D-2
  77.       IF(A-1.D0)6,6,5
  78.     5 HH=HH*A
  79.     6 IF(Z-1.D0)8,8,7
  80.     7 HH=HH/Z
  81.     8 IF(HH-C)10,10,9
  82.     9 HH=C
  83. C
  84. C     INITIALIZE DIFFERENTIATION LOOP
  85.    10 Z=(FCT(X+HH)-FCT(X-HH))/(HH+HH)
  86.       J=5
  87.       JJ=J-1
  88.       AUX(J)=Z
  89.       DH=HH/DFLOAT(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)-FCT(X-HH))/(HH+HH)
  97. C
  98. C     INITIALIZE EXTRAPOLATION LOOP
  99.       D2=1.E38
  100.       B=0.D0
  101.       A=1.D0/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*(2.D0+B))
  108.       AUX(I+1)=AUX(I)+HH
  109. C
  110. C        TEST ON OSCILLATING INCREMENTS
  111.       D2=DABS(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.       END OF EXTRAPOLATION LOOP
  130.