'***************************************************** '* Program to demonstrate the real domain * '* Zbrent subroutine * '* ------------------------------------------------- * '* Reference: BORLAND MATHEMATICAL LIBRARY * '* * '* Basic version by J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '* ------------------------------------------------- * '* Example: Find a real root of f(x)=(x+1)^5 * '* * '* SAMPLE RUN: * '* * '* Input interval (X1,X2): * '* * '* X1 = -2 * '* X2 = 0 * '* * '* Convergence criterion: 1e-10 * '* Maximum number of iterations: 10 * '* * '* The estimated root is: * '* * '* X = -1 * '* * '* The associated Y value is Y = 0 * '* * '* The number of iterations was: 2 * '* The error code is: 0 * '* * '***************************************************** DEFINT I-N DEFDBL A-H, O-Z CLS PRINT PRINT " Input interval (X1,X2):" PRINT INPUT " X1 = ", x1 INPUT " X2 = ", x2 PRINT INPUT " Convergence criterion: ", e INPUT " Maximum number of iterations: ", maxiter GOSUB 2000 'Call BrentRoots subroutine PRINT PRINT " The estimated root is:" PRINT PRINT " X = "; result PRINT PRINT " The associated Y value is Y = "; yr PRINT PRINT " The number of iterations was: "; niter PRINT " The error code is: "; ierror PRINT END 1000 'Test function for Brent method AFunction = (x + 1)^5 RETURN 1200 'TRUE if FA*FB negative IF FA * FB < 0 THEN IRootBracketed = 1 ELSE IRootBracketed = 0 END IF RETURN 1400 'returns the minimum of two real numbers xx1 and xx2 IF xx1 <= xx2 THEN xminimum = xx1 ELSE xminimum = xx2 END IF RETURN '***************************************************** '* Brent Method Function * '* ------------------------------------------------- * '* The purpose is to find a real root of a real * '* function f(x) using Brent method. * '* * '* INPUTS: x1,x2 : interval of root * '* Tolerance : desired accuracy for root * '* maxIter : maximum number of iterations * '* * '* OUTPUTS: The function returns the root value * '* ValueAtRoot : value of f(root) * '* niter : number of done iterations * '* ierror : =0, all OK * '* : =1, no root found in interval * '* : =2, no more iterations ! * '***************************************************** 2000 FPP = 1E-11 xnearzero = 1E-20 i = 0: ierror = 0 'Get values of test function for x=x1 and x=x2 x = x1: GOSUB 1000: FA = AFunction x = x2: GOSUB 1000: FB = AFunction 'Test if a root exists in the interval (x1,x2) GOSUB 1200 'f(x1)*f(x2) not negative? IF IRootBracketed = 0 THEN ierror = 1: RETURN END IF FC = FB : AA = x1 : BB = x2 2100 'main loop If FB*FC > 0 THEN CC=AA : FC=FA : DD = BB - AA : EE = DD END IF IF ABS(FC) < ABS(FB) THEN AA = BB: BB = CC: CC = AA FA = FB: FB = FC: FC = FA END IF Tol1 = 2# * FPP * ABS(BB) + .5# * e xm = .5# * (CC - BB) IF (ABS(xm) <= Tol1) OR (ABS(FA) < xnearzero) THEN result = BB x = result: GOSUB 1000: yr = AFunction GOTO 2500 END IF IF (ABS(EE) >= Tol1) AND (ABS(FA) > ABS(FB)) THEN SS = FB / FA IF ABS(AA - CC) < xnearzero THEN PP = 2# * xm * SS QQ = 1# - SS ELSE QQ = FA / FC RR = FB / FC PP = SS * (2# * xm * QQ * (QQ - RR) - (BB - AA) * (RR - 1#)) QQ = (QQ - 1#) * (RR - 1#) * (SS - 1#) END IF IF PP > xnearzero THEN QQ = -QQ PP = ABS(PP) xx1 = 3# * xm * QQ - ABS(Tol1 * QQ): xx2 = ABS(EE * QQ): GOSUB 1400 IF 2# * PP < xminimum THEN EE = DD: DD = PP / QQ ELSE DD = xm: EE = DD END IF ELSE DD = xm EE = DD END IF AA = BB FA = FB IF ABS(DD) > Tol1 THEN BB = BB + DD ELSE IF xm > 0 THEN BB = BB + ABS(Tol1) ELSE BB = BB - ABS(Tol1) END IF END IF x = BB: GOSUB 1000: FB = AFunction i = i+1 IF i >= maxiter THEN ierror = 2: GOTO 2500 END IF 'Continue GOTO 2100 2500 'A real root has been found niter = i RETURN 'End of file zbrent.bas