home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / fortran / library / ssp / rootnleq / drtmi.for next >
Encoding:
Text File  |  1987-12-06  |  5.0 KB  |  166 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE DRTMI
  5. C
  6. C        PURPOSE
  7. C           TO SOLVE GENERAL NONLINEAR EQUATIONS OF THE FORM FCT(X)=0
  8. C           BY MEANS OF MUELLER-S ITERATION METHOD.
  9. C
  10. C        USAGE
  11. C           CALL DRTMI (X,F,FCT,XLI,XRI,EPS,IEND,IER)
  12. C           PARAMETER FCT REQUIRES AN EXTERNAL STATEMENT.
  13. C
  14. C        DESCRIPTION OF PARAMETERS
  15. C           X      - DOUBLE PRECISION RESULTANT ROOT OF EQUATION
  16. C                    FCT(X)=0.
  17. C           F      - DOUBLE PRECISION RESULTANT FUNCTION VALUE
  18. C                    AT ROOT X.
  19. C           FCT    - NAME OF THE EXTERNAL DOUBLE PRECISION FUNCTION
  20. C                    SUBPROGRAM USED.
  21. C           XLI    - DOUBLE PRECISION INPUT VALUE WHICH SPECIFIES THE
  22. C                    INITIAL LEFT BOUND OF THE ROOT X.
  23. C           XRI    - DOUBLE PRECISION INPUT VALUE WHICH SPECIFIES THE
  24. C                    INITIAL RIGHT BOUND OF THE ROOT X.
  25. C           EPS    - SINGLE PRECISION INPUT VALUE WHICH SPECIFIES THE
  26. C                    UPPER BOUND OF THE ERROR OF RESULT X.
  27. C           IEND   - MAXIMUM NUMBER OF ITERATION STEPS SPECIFIED.
  28. C           IER    - RESULTANT ERROR PARAMETER CODED AS FOLLOWS
  29. C                     IER=0 - NO ERROR,
  30. C                     IER=1 - NO CONVERGENCE AFTER IEND ITERATION STEPS
  31. C                             FOLLOWED BY IEND SUCCESSIVE STEPS OF
  32. C                             BISECTION,
  33. C                     IER=2 - BASIC ASSUMPTION FCT(XLI)*FCT(XRI) LESS
  34. C                             THAN OR EQUAL TO ZERO IS NOT SATISFIED.
  35. C
  36. C        REMARKS
  37. C           THE PROCEDURE ASSUMES THAT FUNCTION VALUES AT INITIAL
  38. C           BOUNDS XLI AND XRI HAVE NOT THE SAME SIGN. IF THIS BASIC
  39. C           ASSUMPTION IS NOT SATISFIED BY INPUT VALUES XLI AND XRI, THE
  40. C           PROCEDURE IS BYPASSED AND GIVES THE ERROR MESSAGE IER=2.
  41. C
  42. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  43. C           THE EXTERNAL DOUBLE PRECISION FUNCTION SUBPROGRAM FCT(X)
  44. C           MUST BE FURNISHED BY THE USER.
  45. C
  46. C        METHOD
  47. C           SOLUTION OF EQUATION FCT(X)=0 IS DONE BY MEANS OF MUELLER-S
  48. C           ITERATION METHOD OF SUCCESSIVE BISECTIONS AND INVERSE
  49. C           PARABOLIC INTERPOLATION, WHICH STARTS AT THE INITIAL BOUNDS
  50. C           XLI AND XRI. CONVERGENCE IS QUADRATIC IF THE DERIVATIVE OF
  51. C           FCT(X) AT ROOT X IS NOT EQUAL TO ZERO. ONE ITERATION STEP
  52. C           REQUIRES TWO EVALUATIONS OF FCT(X). FOR TEST ON SATISFACTORY
  53. C           ACCURACY SEE FORMULAE (3,4) OF MATHEMATICAL DESCRIPTION.
  54. C           FOR REFERENCE, SEE G. K. KRISTIANSEN, ZERO OF ARBITRARY
  55. C           FUNCTION, BIT, VOL. 3 (1963), PP.205-206.
  56. C
  57. C     ..................................................................
  58. C
  59.       SUBROUTINE DRTMI(X,F,FCT,XLI,XRI,EPS,IEND,IER)
  60. C
  61. C
  62.       DOUBLE PRECISION X,F,FCT,XLI,XRI,XL,XR,FL,FR,TOL,TOLF,A,DX,XM,FM
  63. C
  64. C     PREPARE ITERATION
  65.       IER=0
  66.       XL=XLI
  67.       XR=XRI
  68.       X=XL
  69.       TOL=X
  70.       F=FCT(TOL)
  71.       IF(F)1,16,1
  72.     1 FL=F
  73.       X=XR
  74.       TOL=X
  75.       F=FCT(TOL)
  76.       IF(F)2,16,2
  77.     2 FR=F
  78.       IF(DSIGN(1.D0,FL)+DSIGN(1.D0,FR))25,3,25
  79. C
  80. C     BASIC ASSUMPTION FL*FR LESS THAN 0 IS SATISFIED.
  81. C     GENERATE TOLERANCE FOR FUNCTION VALUES.
  82.     3 I=0
  83.       TOLF=100.*EPS
  84. C
  85. C
  86. C     START ITERATION LOOP
  87.     4 I=I+1
  88. C
  89. C     START BISECTION LOOP
  90.       DO 13 K=1,IEND
  91.       X=.5D0*(XL+XR)
  92.       TOL=X
  93.       F=FCT(TOL)
  94.       IF(F)5,16,5
  95.     5 IF(DSIGN(1.D0,F)+DSIGN(1.D0,FR))7,6,7
  96. C
  97. C     INTERCHANGE XL AND XR IN ORDER TO GET THE SAME SIGN IN F AND FR
  98.     6 TOL=XL
  99.       XL=XR
  100.       XR=TOL
  101.       TOL=FL
  102.       FL=FR
  103.       FR=TOL
  104.     7 TOL=F-FL
  105.       A=F*TOL
  106.       A=A+A
  107.       IF(A-FR*(FR-FL))8,9,9
  108.     8 IF(I-IEND)17,17,9
  109.     9 XR=X
  110.       FR=F
  111. C
  112. C     TEST ON SATISFACTORY ACCURACY IN BISECTION LOOP
  113.       TOL=EPS
  114.       A=DABS(XR)
  115.       IF(A-1.D0)11,11,10
  116.    10 TOL=TOL*A
  117.    11 IF(DABS(XR-XL)-TOL)12,12,13
  118.    12 IF(DABS(FR-FL)-TOLF)14,14,13
  119.    13 CONTINUE
  120. C     END OF BISECTION LOOP
  121. C
  122. C     NO CONVERGENCE AFTER IEND ITERATION STEPS FOLLOWED BY IEND
  123. C     SUCCESSIVE STEPS OF BISECTION OR STEADILY INCREASING FUNCTION
  124. C     VALUES AT RIGHT BOUNDS. ERROR RETURN.
  125.       IER=1
  126.    14 IF(DABS(FR)-DABS(FL))16,16,15
  127.    15 X=XL
  128.       F=FL
  129.    16 RETURN
  130. C
  131. C     COMPUTATION OF ITERATED X-VALUE BY INVERSE PARABOLIC INTERPOLATION
  132.    17 A=FR-F
  133.       DX=(X-XL)*FL*(1.D0+F*(A-TOL)/(A*(FR-FL)))/TOL
  134.       XM=X
  135.       FM=F
  136.       X=XL-DX
  137.       TOL=X
  138.       F=FCT(TOL)
  139.       IF(F)18,16,18
  140. C
  141. C     TEST ON SATISFACTORY ACCURACY IN ITERATION LOOP
  142.    18 TOL=EPS
  143.       A=DABS(X)
  144.       IF(A-1.D0)20,20,19
  145.    19 TOL=TOL*A
  146.    20 IF(DABS(DX)-TOL)21,21,22
  147.    21 IF(DABS(F)-TOLF)16,16,22
  148. C
  149. C     PREPARATION OF NEXT BISECTION LOOP
  150.    22 IF(DSIGN(1.D0,F)+DSIGN(1.D0,FL))24,23,24
  151.    23 XR=X
  152.       FR=F
  153.       GO TO 4
  154.    24 XL=X
  155.       FL=F
  156.       XR=XM
  157.       FR=FM
  158.       GO TO 4
  159. C     END OF ITERATION LOOP
  160. C
  161. C
  162. C     ERROR RETURN IN CASE OF WRONG INPUT DATA
  163.    25 IER=2
  164.       RETURN
  165.       END
  166.