'********************************************************************* '* Find a real root of a real function F(x) by the Zeroin method * '* ----------------------------------------------------------------- * '* SAMPLE RUN: * '* (Find a root of function 5x - exp(x) between 0 and 1) * '* * '* The output file tzeroin.lst contains: * '* * '* ---------------------------------------------------- * '* Zeroin method for real valued, nonlinear functions * '* ---------------------------------------------------- * '* Find a real root of function: * '* F(x) = 5x - Exp(x) * '* Starting value x0 = 0 * '* Return code = 0 * '* Root = .259171101819074 * '* Function value = 1.113367210925365D-15 * '* Function calls = 8 * '* Absolute error = 1.20000004670464D-14 * '* Relative error = 1.20000004670464D-13 * '* ---------------------------------------------------- * '* * '* ----------------------------------------------------------------- * '* Ref.: "Numerical Algorithms with C By G. Engeln-Mueller and * '* F. Uhlig, Springer-Verlag, 1996" [BIBLI 11]. * '* * '* Basic version by J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '********************************************************************* DEFDBL A-H, O-Z DEFINT I-N FOUR = 4# HALF = .5# ONE = 1# THREE = 3# TWO = 2# ZERO = 0# XMACHEPS = 1.2E-16 'machine smallest real number ' the function Fkt(x,y) is defined in subroutine 100 abserr = 100# * XMACHEPS relerr = 1000# * XMACHEPS ' open output file OPEN "tzeroin.lst" FOR OUTPUT AS #1 CLS PRINT #1, "----------------------------------------------------" PRINT #1, " Zeroin method for real valued, nonlinear functions " PRINT #1, "----------------------------------------------------" ifmax = 100: a = ZERO: b = ONE PRINT #1, " Find a real root of function:" PRINT #1, " F(x) = 5x - Exp(x)" PRINT #1, " Starting value x0 = ", a GOSUB 1000 'call zeroin(abserr, relerr, ifmax, a, b, fb, ifanz, irc) PRINT #1, " Return code = "; irc PRINT #1, " Root = "; b PRINT #1, " Function value = "; fb PRINT #1, " Function calls = "; ifanz PRINT #1, " Absolute error = "; abserr PRINT #1, " Relative error = "; relerr PRINT #1, "----------------------------------------------------" CLOSE #1 PRINT PRINT " Results in tzeroin.lst." PRINT END 'test function F(x) 100 'Subroutine Fkt(x,y) y = 5# * x - EXP(x) RETURN 'Subroutine zeroin Find roots with the Zeroin method ' ( ' abserr, absolute error bound .............. ' relerr, relative error bound .............. ' ifmax, maximal number of calls for fkt() . ' a, [a,b] = inclusion interval ........ ' b, right endpoint or zero ............ ' fb, Function value at the root b ...... ' ifanz, number of actual function calls ... ' irc error code ........................ ' ) '************************************************************************ '* Given a real valued function f on an interval [a,b] with * '* f(a) * f(b) < 0, we compute a root of f in [a,b]. * '* The Zeroin method combines bisection and secant methods with inverse * '* quadratic interpolation. * '* * '* Input parameters: * '* ================= * '* abserr\ Error bounds with abserr >= 0 and relerr >= 0. Their * '* relerr/ sum must exceed zero. For break-off we use the test * '* |xm| <= 0.5 * (abserr + |b| * relerr), * '* where xm denotes half the length of the final inclusion * '* interval. For relerr = 0 we test only for absolute error, * '* while for abserr = 0, we test the relative error only. * '* abserr and relerr are used as put in only when both exceed * '* four times the machine constant. Otherwise we set them to * '* this value. * '* ifmax upper limit of calls of function * '* a,b end points of the interval, that includes a root * '* * '* Output parameters: * '* ================= * '* abserr\ actual error bounds used * '* relerr/ * '* b approximate root * '* fb value of f at the root b * '* ifanz number of calls of fkt() * '* * '* Error code irc: * '* ============== * '* = -2: abserr or relerr is negative, or both are zero, or ifmax < 1 * '* = -1: The necessary assumption f(a) * f(b) < 0 is violated. * '* = 0: Desired accuracy has been reache : * '* |xm| <= 0.5 * (abserr + |b| * relerr). * '* = 1: b on output is a root with fkt(b) = 0. * '* = 2: Either a or b is a root at input. * '* = 3: After fmax calls of Fkt the algorithm has not found a root * '* = 4: Error when opening intermediate result file * '* * '************************************************************************ 'Labels: 1010, 1050 'double fa, & " Function value fkt(a) ' fc, & " Function value fkt(c) ' eps, & " minimal value for error bounds abserr and relerr ' tol1, & " auxiliary variable for mixed error bound ' xm, & " half of the current interval length ' c, & " value inside [a,b] ' d, & " Distance to the nearest approximate root ' e, & " previous value of d ' p, & ' q, & ' r, & ' s, & ' tmp auxiliary variable to check inclusion ' f(a) * f(b) < 0 'integer ifehler error code of this function ' ------------------ initialize variables ----------------------------- 1000 c = ZERO: d = ZERO: e = ZERO: irc = 0: ifehler = 0: eps = FOUR * XMACHEPS x = a: GOSUB 100 'call Fkt(a,fa) evaluate fkt() at the end fa = y 'points a and b x = b: GOSUB 100 'call Fkt(b,fb) fb = y ifanz = 2 ' ----------- check f(a) * f(b) < 0 -------------------------------- tmp = fa * fb IF tmp > ZERO THEN irc = -1 ELSEIF tmp = ZERO THEN irc = 2 END IF IF irc <> 0 THEN RETURN ' ----------- check usability of given error bounds ------------------- IF abserr < ZERO OR relerr < ZERO OR abserr + relerr <= ZERO OR ifmax < 1 THEN rc = -2 RETURN END IF IF relerr = ZERO AND abserr < eps THEN abserr = eps ELSEIF abserr = ZERO AND relerr < eps THEN relerr = eps ELSE IF abserr < eps THEN abserr = eps IF relerr < eps THEN relerr = eps END IF ' open file for detailed results OPEN "tzeroin.log" FOR OUTPUT AS #2 fc = fb ' start iteration 1010 IF fb * (fc / ABS(fc)) > ZERO THEN ' no inclusion of a root c = a ' between b and c ? fc = fa ' alter c so that b and c e = b - a ' include the root of f d = e END IF IF ABS(fc) < ABS(fb) THEN ' If fc has the smaller modulus a = b ' interchange interval end b = c ' points c = a fa = fb fb = fc fc = fa END IF PRINT #2, " a ="; a; " b ="; b; " c ="; c PRINT #2, " fa="; fa; " fb="; fb; " fc="; fc tol1 = HALF * (abserr + relerr * ABS(b)) xm = HALF * (c - b) PRINT #2, " xm="; xm; " tol1="; tol1 IF ABS(xm) <= tol1 THEN ' reached desired accuracy ? ifehler = 0 RETURN ' normal exit END IF IF fb = ZERO THEN ' Is the best approximate root b ifehler = 1 ' a precise root of f ? RETURN ' other normal exit, b is a root END IF r = ZERO IF ABS(e) < tol1 OR ABS(fa) <= ABS(fb) THEN e = xm d = e PRINT #2, " Bisection" ELSE IF a <> c THEN ' if a is not equal to c q = fa / fc ' With a, b and c we have 3 points for r = fb / fc ' an inverse quadratic interpolation s = fb / fa p = s * (TWO * xm * q * (q - r) - (b - a) * (r - ONE)) q = (q - ONE) * (r - ONE) * (s - ONE) ELSE ' If a equals c : s = fb / fa ' Use the secant method or linear p = TWO * xm * s ' interpolation q = ONE - s END IF IF p > ZERO THEN ' alter the sign of p/q for the q = -q ' subsequent division ELSE p = ABS(p) END IF IF TWO * p >= THREE * xm * q - ABS(tol1 * q) OR p >= ABS(HALF * e * q) THEN e = xm: d = e PRINT #2, " Bisection" ELSE e = d ' Compute the quotient p/q for both iterations d = p / q ' which will be used to modify b IF r = ZERO THEN PRINT #2, " Secant method" ELSE PRINT #2, " Inverse quadratic interpolation" END IF END IF END IF a = b ' store the best approximation b and its fa = fb ' function value fb in a and fa IF ABS(d) > tol1 THEN ' d large enough? b = b + d ' add d to b PRINT #2, " Difference d from new b: "; d ELSE ' d too small? IF xm > ZERO THEN tmp = ABS(tol1) ELSE tmp = -ABS(tol1) END IF b = b + tmp ' add tol1 to b IF xm < ZERO THEN PRINT #2, " Subtract error bound: d = "; -tol1 ELSE PRINT #2, " Add error bound: d = "; tol1 END IF END IF x = b GOSUB 100 'call Fkt(b,fb) ' compute function value at b fb = y ifanz = ifanz + 1 ' up evaluation counter PRINT #2, " b = "; b; " fb= "; fb; " Number of functional evaluations = "; ifanz IF ifanz > ifmax THEN ' too many function evaluations? ifehler = 3 GOTO 1050 END IF GOTO 1010 ' end iteration 1050 CLOSE #2 irc = ifehler RETURN 'end of file tzeroin.bas