SUBROUTINE MRF(FUNC,XLO,XHI,X0,ERROR,TOL,MAX,IERR) C C This subroutine uses modified regula falsi to find a root of C the user supplied function. The root must satisfy both an C error bound on X0 and a function tolerance on FUNC( X0 ). C Extrapolation is performed when necessary. FCTR: MRF C C FUNC Function name, must be declared EXTERNAL in calling C program. Assumed to be a real function of a single C real variable. C XLO Lower bound on root. C XHI Upper bound on root. C X0 Best approximation of root. C ERROR Allowable error on X0. C TOL Allowable error on FUNC( X0 ). C MAX Maximum number of iterations to perform. C IERR Return code, where: C = 0, Normal return. C = -1, Did not satisfy error requirement on X0. C = -2, Did not satisfy tolerance requirement on C FUNC( X0 ). C = -4, Did not find root in MAX iterations. C = -8, Extrapolation, X0 < XLO. C = -16, Extrapolation, X0 > XHI. C Note: IERR may be the sum of any of the above. C C C - INITIALIZE FOR ROOT SEARCH. F = FUNC(XLO) G = FUNC(XHI) A = XLO B = XHI W = A FA = F FW = F C C - GET SET TO DO NO MORE THAN MAX ITERATIONS. DO 1000 I=1,MAX C C - CALCULATE A NEW ESTIMATE FOR X0. WP = (G*A-F*B)/(G-F) FWP = FUNC(WP) C C - IS NEW ROOT ON SAME SIDE AS OLD ROOT? 400 IF( FA*FWP .GT. 0.0 ) GO TO 500 B = WP G = FWP IF( FW*FWP .GT. 0.0 ) F = F*0.5 GO TO 600 500 A = WP F = FWP IF( FW*FWP .GT. 0.0 ) G = G*0.5 C C - SATISFIES BOTH CRITERIA FOR A ROOT? 600 IF( ABS(B-A) .GT. ERROR ) GO TO 1000 IF( ABS(FWP) .LE. TOL ) GO TO 1100 1000 FW = FWP IERR = -4 IF( ABS(B-A) .GT. ERROR ) IERR = IERR - 1 IF( ABS(FWP) .GT. TOL ) IERR = IERR - 2 1100 X0 = (A+B)*0.5 IF( X0 .LT. XLO ) IERR = IERR - 8 IF( X0 .GT. XHI ) IERR = IERR - 16 1200 RETURN END