/**************************************************** * Program to demonstrate the real domain * * Zbrent subroutine * * ------------------------------------------------- * * Reference: BORLAND MATHEMATICAL LIBRARY * * * * C++ 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.000000 * * * * The associated Y value is Y = 0.000000 * * * * The number of iterations was: 2 * * The error code is: 0 * * * ****************************************************/ #include #include #define FALSE 0 #define TRUE 1 double e,r,x1,x2,yr; int err,maxiter,k; // Test function for Brent method double AFunction(double x) { return (x+1)*(x+1)*(x+1)*(x+1)*(x+1); } // TRUE if x1*x2 negative int RootBracketed(double x1,double x2) { int result; if ((x1 > 0 && x2 > 0) || (x1 < 0 && x2 < 0)) result = FALSE; else result = TRUE; return result; } // returns the minimum of two real numbers double Minimum(double x1,double x2) { double result; if (x1 < x2) result = x1; else result = x2; return result; } /**************************************************** * Brent Method Function * * ------------------------------------------------- * * The purpose is to find a real root of a real * * function f(x) using Brent's 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 * * error : =0, all OK * * : =1, no root found in interval * * : =2, no more iterations ! * ****************************************************/ double BrentRoots( double x1, double x2, double Tolerance, int maxIterations, double *valueAtRoot, int *niter, int *error ) { double FPP = 1e-11, nearzero = 1e-20; double result, AA, BB, CC, DD, EE, FA, FB, FC, Tol1, PP, QQ, RR, SS, xm; int i, done; i = 0; done = FALSE; error = 0; AA = x1; BB = x2; FA = AFunction(AA); FB = AFunction(BB); if (!(RootBracketed(FA,FB))) *error = 1; else { FC = FB; do { if (!(RootBracketed(FC,FB))) { CC = AA; FC = FA; DD = BB - AA; EE = DD; } if (fabs(FC) < fabs(FB)) { AA = BB; BB = CC; CC = AA; FA = FB; FB = FC; FC = FA; } Tol1 = 2.0 * FPP * fabs(BB) + 0.5 * Tolerance; xm = 0.5 * (CC-BB); if ((fabs(xm) <= Tol1) || (fabs(FA) < nearzero)) { result = BB; done = TRUE; *valueAtRoot = AFunction(result); } // A root has been found else { if ((fabs(EE) >= Tol1) && (fabs(FA) > fabs(FB))) { SS = FB/ FA; if (fabs(AA - CC) < nearzero) { PP = 2.0 * xm * SS; QQ = 1.0 - SS; } else { QQ = FA/FC; RR = FB /FC; PP = SS * (2.0 * xm * QQ * (QQ - RR) - (BB-AA) * (RR - 1.0)); QQ = (QQ - 1.0) * (RR - 1.0) * (SS - 1.0); } if (PP > nearzero) QQ = -QQ; PP = fabs(PP); if ((2.0 * PP) < Minimum(3.0*xm *QQ-fabs(Tol1 * QQ), fabs(EE * QQ))) { EE = DD; DD = PP/QQ; } else { DD = xm; EE = DD; } } else { DD = xm; EE = DD; } AA = BB; FA = FB; if (fabs(DD) > Tol1) BB = BB + DD; else { if (xm > 0) BB = BB + fabs(Tol1); else BB = BB - fabs(Tol1); } FB = AFunction(BB); i++; } } while ((!done) && (i < maxIterations)); if (i >= maxIterations) *error = 2; } *niter = i; return result; } // BrentRoots() void main() { printf("\n Input interval (X1,X2):\n\n"); printf(" X1 = "); scanf("%lf",&x1); printf(" X2 = "); scanf("%lf",&x2); printf("\n Convergence criterion: "); scanf("%lf",&e); printf("\n Maximum number of iterations: "); scanf("%d",&maxiter); r = BrentRoots(x1,x2,e,maxiter,&yr,&k,&err); printf("\n\n The estimated root is:\n\n"); printf(" X = %f\n\n", r); printf(" The associated Y value is Y = %f\n\n", yr); printf(" The number of iterations was: %d\n",k); printf(" The error code is: %d\n\n",err); } // End of file zbrent.cpp