/**************************************************************** * Program to demonstrate Bisection & Quartile subroutines * * ------------------------------------------------------------- * * References: * * * * Bisection: BASIC Scientific Subroutines, Vol. II * * By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981. * * * * Quartile: After an algorithm provided By Namir Shammas. * * * * C++ Release By J-P Moreau, Paris. * * (www.jpmoreau.fr) * * ------------------------------------------------------------- * * Example: Find a real root of f(x)=(x+1)^5 * * * * Sample run: * * * * What is the initial range (X0,X1): * * X0 = -5 * * X1 = 0 * * Convergence criterion: 1e-6 * * * * * * Bisection: The calculated zero is X = -9.999996e-001 * * * * The associated Y value is Y = 5.850012e-033 * * * * The number of steps was: 23 * * * * * * Quartile: The calculated zero is X = -1.000000e+000 * * * * The associated Y value is Y = -5.712903e-036 * * * * The number of steps was: 12 * * * ***************************************************************** Note: The quartile method is slightly faster and more accurate than the classical Bisection Method. ---------------------------------------------------------------*/ #include #include double a,b,root, e,x,x0,x1; int m; /**********************************************/ double Y(double x) { return (x+1)*(x+1)*(x+1)*(x+1)*(x+1); } /**********************************************/ /********************************************** * Bisection method subroutine * * ------------------------------------------- * * This routine iteratively seeks the zero of * * function Y(x) using the method of interval * * halving until the interval is less than e * * in width. It is assumed that the function * * Y(x) is available from a function routine. * * ------------------------------------------- * * Input values: range (x0,x1), and e. * * Output values: root x, Y(x) and number of * * steps m. * **********************************************/ void Bisect(void) { //Labels: e10 double y0, yy; m=0; if (Y(x0) * Y(x1) > 0.0) { printf(" Bisection: No root found in given interval.\n"); return; } else { e10: y0=Y(x0); x=(x0+x1)/2; yy=Y(x); m++; if (yy*y0==0.0) return; if (yy*y0<0.0) x1=x; if (yy*y0>0.0) x0=x; if (fabs(x1-x0)>e) goto e10; } } /*********************************************** * Quartile method subroutine * * -------------------------------------------- * * This routine iteratively seeks the zero of * * function Y(x) using the method of interval * * quarting until the interval is less than tol * * in width. It is assumed that the function * * Y(x) is available from a function routine. * * -------------------------------------------- * * Input values: range (a, b), and tol. * * Output values: found root and number of used * * steps. * ***********************************************/ void Quartile(double a,double b,double tol, double *root, int *step) { double m, co; co=0.25; //can be 0.3 *step=0; if (Y(a) * Y(b) > 0.0) { printf(" Quartile: No root found in given interval.\n"); return; } else { while (fabs(a-b) > tol) { if (fabs(Y(a)) < fabs(Y(b))) m = a + co*(b-a); else m = a + (1.0-co)*(b-a); if (Y(m)*Y(a) > 0.0) a= m; else b=m; (*step)=(*step)+1; } *root = (a+b)/2.0; } } void main() { printf("\n What is the initial range (X0,X1):\n\n"); printf(" X0 = "); scanf("%lf", &x0); printf(" X1 = "); scanf("%lf", &x1); printf("\n Convergence criterion: "); scanf("%lf", &e); printf("\n"); a=x0; b=x1; Bisect(); // Call Bisection routine if (m!=0) { printf("\n Bisection: The calculated zero is X = %e\n\n", x); printf(" The associated Y value is Y = %e\n\n", Y(x)); printf(" The number of steps was: %d\n\n", m); } Quartile(a, b, e, &root, &m); //Call Quartile routine if (m!=0) { printf("\n Quartile: The calculated zero is X = %e\n\n", root); printf(" The associated Y value is Y = %e\n\n", Y(root)); printf(" The number of steps was: %d\n\n", m); } } //End of file tquart.cpp