/**************************************************** * Program to demonstrate Mueller's method * * ------------------------------------------------- * * Reference: BASIC Scientific Subroutines, Vol. II * * By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1]. * * * * C++ version by J-P Moreau, Paris. * * (www.jpmoreau.fr) * * ------------------------------------------------- * * Example: Find a real root of f(x)=(x+1)^5 * * * * Sample run: * * * * Input the initial guess: * * X0 = 0 * * What is the bond of this guess: 3 * * Error criterion: 1e-6 * * Maximum number of iterations: 100 * * * * The calculated zero is X = -1.000476 * * The associated value is Y = 0.000000 * * The number of steps was: 47 * * * ****************************************************/ #include #include // global variables int k,n; double d,e,x,x0; //*************************************************** // Function subroutine double Y(double x) { return (1+5*x+10*x*x+10*x*x*x+5*x*x*x*x+x*x*x*x*x); } //*************************************************** /********************************************** * Mueller's method subroutine * * ------------------------------------------- * * This routine iteratively seeks the root of * * a function Y(x) by fitting a parabola to * * three points and calculating the nearest * * root as described in Becket and Hurt, nume- * * rical calculations and algorithms. * * writeS: * * The routine requires an initial guess, x0, * * a bound on the error on this guess, d and a * * convergence criteria, e. Also required is a * * limit on the number of iterations, n. * * OUTPUTS: * * The routine returns the value of the root * * found, x and the number of iterations per- * * formed, k. * **********************************************/ void Mueller() { //Labels: e100,e200 double a1,b,c1,d1,e1,e2,e3,x1,x2,x3,xl,xl1; // Set up the three evaluation points k=1; x3=x0; x1=x3-d; x2=x3+d; // Calculate Mueller parameters and guard against divide by zero if (x2-x1==0) x2=x2*1.0000001; e100: if (x2-x1==0) x2=x2+0.0000001; xl1=(x3-x2)/(x2-x1); d1=(x3-x1)/(x2-x1); if (k>1) goto e200; // Get values of function e1=Y(x1); e2=Y(x2); e200: e3=Y(x3); a1=xl1*xl1*e1-d1*d1*e2+(xl1+d1)*e3; c1=xl1*(xl1*e1-d1*e2+e3); b=a1*a1-4.0*d1*c1*e3; // Test for complex root, meaning the parabola is inverted if (b<0) b=0.0; // Choose closest root if (a1<0) a1=a1-sqrt(b); if (a1>0) a1=a1+sqrt(b); // Guard against divide by zero if (a1==0) a1=1e-7; xl=-2.0*d1*e3/a1; // Calculate next estimate x=x3+xl*(x3-x2); // Test for convergence if (fabs(x-x3)=n) return; // Otherwise, make another pass k++; // Some saving x1=x2; x2=x3; x3=x; e1=e2; e2=e3; goto e100; } // Mueller() void main() { printf("\n Input the initial guess:\n"); printf("\n X0 = "); scanf("%lf",&x0); printf("\n What is the bond of this guess: "); scanf("%lf",&d); printf("\n Error criterion: "); scanf("%lf",&e); printf("\n Maximum number of iterations: "); scanf("%d",&n); Mueller(); // Call Mueller's routine printf("\n\n The calculated zero is X = %f\n", x); printf("\n The associated Y value is Y = %f\n", Y(x)); printf("\n The number of steps was: %d\n\n", k); } // End of file Mueller.cpp