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