/****************************************************** * Program to demonstrate two dimensional * * 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,y)=(x+1)^5*(y-1)^5 * * * * Sample run: * * * * Input the initial guesses and their bonds: * * X0 = 0 * * Bond on X0 = 3 * * Y0 = 0 * * Bond on Y0 = 3 * * Error criterion: 1e-6 * * Maximum number of iterations: 100 * * * * The calculated zero is X= -0.999998 Y= 0.999998 * * The associated W value is W = -0.000000 * * The number of steps was: 43 * * * ******************************************************/ #include #include // global variables int k,n; double b1,b2,e,x,y,x0,yy0; // shared by Utility() and Mueller2D() double a1,b,c1,d1,e1,e2,e3,xl,xl1; //********************************************** // Function subroutine double W(double x, double y) { double temp; temp = (x+1)*(x+1)*(x+1)*(x+1)*(x+1); return (temp * (y-1)*(y-1)*(y-1)*(y-1)*(y-1)); } //********************************************** // Utility function used by Mueller2D() void Utility() { 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; // Choose closest root if (a1<0) a1=a1-sqrt(b); if (a1>0) a1=a1+sqrt(b); // Guard against a divide by zero if (fabs(a1)+fabs(b)==0) a1=4.0*d1*e3; // Calculate a relative distance of next gues // and guard against a divide by zero if (a1==0) a1=1e-7; xl=-2.0*d1*e3/a1; } /************************************************ * Two Dimensional Mueller's method subroutine * * --------------------------------------------- * * This routine iteratively seeks the root of a * * function W(x,y) by fitting a parabola to * * three points and calculating the nearest root * * as described in Becket and Hurt, numerical * * calculations and algorithms. * * INPUTS: * * x0, y0 - The initial guess * * b1, b2 - A bound on the error in this guess * * e - The convergence criteria * * n - The maximum number of iterations * * OUTPUTS: * * x,y - The value of the root found * * k - The number of iterations performed. * ************************************************/ void Mueller2D() { //Label: e100 double x1,x2,x3,y1,y2,y3; // Set up the three evaluation points k=1; e100: x3=x0; x1=x3-b1; x2=x3+b1; // Calculate Mueller parameters and guard against divide by zero if (x2-x1==0) x2=x2*1.0000001; if (x2-x1==0) x2=x2+0.0000001; xl1=(x3-x2)/(x2-x1); d1=(x3-x1)/(x2-x1); // Get values of function e1=W(x1,yy0); e2=W(x2,yy0); e3=W(x3,yy0); Utility(); // Calculate new x estimate b1=xl*(x3-x2); x=x3+b1; // Test for convergence if (fabs(b1)+fabs(b2)=n) return; yy0=y; k++; // Start another pass goto e100; } //Mueller2D() void main() { printf("\n Input the initial guess:\n\n"); printf(" X0 = "); scanf("%lf",&x0); printf(" Bond on X0 = "); scanf("%lf",&b1); printf("\n Y0 = "); scanf("%lf",&yy0); printf(" Bond on Y0 = "); scanf("%lf",&b2); printf("\n Error criterion: "); scanf("%lf",&e); printf(" Maximum number of iterations: "); scanf("%d",&n); Mueller2D(); // Call 2D Mueller's routine printf("\n\n The calculated zero is (X,Y) = %f %f\n", x, y); printf(" The associated Y value is Y = %f\n", W(x,y)); printf(" The number of steps was: %d\n\n", k); } // End of file Mueller2.cpp