/**************************************************** * Program to demonstrate the complex domain * * Mueller's subroutine * * ------------------------------------------------- * * 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 complex root of zē + 1 = 0 * * * * SAMPLE RUN: * * * * Write the initial guesses and their bounds: * * * * X0 = 0 * * Bond on X0 = 3 * * Y0 = 0 * * Bond on Y0 = 3 * * * * Convergence criterion: 0 * * * * Maximum number of iterations: 100 * * * * The calculated zero is (X,Y) = * * -0.000014 -1.000034 * * * * The associated Z value is (U,V) = * * -0.000067 0.000028 * * * * The number of steps was: 100 * ****************************************************/ #include #include int k,n; double b1,b2,e,u,v,x,y,x0,yy0; //********************************************** // Function subroutine void Z(double x,double y,double *u,double *v) { *u = x*x-y*y+1; *v = 2.0*x*y; } //********************************************** /************************************************ * Mueller's method for complex roots subroutine * * --------------------------------------------- * * This routine uses the parabolic fitting tech- * * nique associated with Mueller's method, but * * does it in the complex domain. * * --------------------------------------------- * * 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 complex root found * * k - The number of iterations performed, * ************************************************/ void ZMueller() { //Labels: e100,e200,e300 double d,d1,d2,x1,x2,x3,y1,y2,y3; double u1,u2,u3,v1,v2,v3,xl1,xl2; double a,a1,a2,b,c1,c2,e1,e2; // Start calculations k=1; x3=x0; y3=yy0; x1=x3-b1; y1=y3-b2; x2=x3+b1; y2=y3+b2; e100: d=(x2-x1)*(x2-x1)+(y2-y1)*(y2-y1); // Avoid divide by zero if (d==0) d=1e-7; xl1=(x3-x2)*(x2-x1)+(y3-y2)*(y2-y1); xl1=xl1/d; xl2=(x2-x1)*(y3-y2)+(x3-x2)*(y2-y1); xl2=xl2/d; d1=(x3-x1)*(x2-x1)+(y3-y1)*(y2-y1); d1=d1/d; d2=(x2-x1)*(y3-y1)+(x3-x1)*(y2-y1); d2=d2/d; // Get function values Z(x1,y1,&u1,&v1); Z(x2,y2,&u2,&v2); Z(x3,y3,&u3,&v3); // Calculate Mueller parameters e1=u1*(xl1*xl1-xl2*xl2)-2.0*v1*xl1*xl2-u2*(d1*d1-d2*d2); e1=e1+2.0*v2*d1*d2+u3*(xl1+d1)-v3*(xl2+d2); e2=2.0*xl1*xl2*u1+v1*(xl1*xl1-xl2*xl2)-2.0*d1*d2*u2-v2*(d1*d1-d2*d2); e2=e2+u3*(xl2+d2)+v3*(xl1+d1); c1=xl1*xl1*u1-xl1*xl2*v1-d1*xl1*u2+xl1*d2*v2+u3*xl1; c1=c1-u1*xl2*xl2-v1*xl1*xl2+u2*xl2*d2+v2*d1*xl2-v3*xl2; c2=xl1*xl2*u1+xl1*xl1*v1-d2*xl1*u2-xl1*d1*v2+v3*xl1; c2=c2+u1*xl1*xl2-v1*xl2*xl2-u2*xl2*d1+v2*d2*xl2+u3*xl2; b1=e1*e1-e2*e2-4.0*(u3*d1*c1-u3*d2*c2-v3*d2*c1-v3*d1*c2); b2=2.0*e1*e2-4.0*(u3*d2*c1+u3*d1*c2+v3*d1*c1-v3*d2*c2); // Guard against divide by zero if (b1==0) b1=1e-7; a=atan(b2/b1); a=a/2.0; b=sqrt(sqrt(b1*b1+b2*b2)); b1=b*cos(a); b2=b*sin(a); a1=(e1+b1)*(e1+b1)+(e2+b2)*(e2+b2); a2=(e1-b1)*(e1-b1)+(e2-b2)*(e2-b2); if (a1>a2) goto e200; a1=e1-b1; a2=e2-b2; goto e300; e200: a1=e1+b1; a2=e2+b2; e300: a=a1*a1+a2*a2; xl1=a1*d1*u3-a1*d2*v3+a2*u3*d2+a2*v3*d1; // Guard against divide by zero if (a==0) a=1e-7; xl1=-2.0*xl1/a; xl2=-d1*u3*a2+d2*v3*a2+a1*u3*d2+a1*v3*d1; xl2=-2.0*xl2/a; // Calculate new estimate x=x3+xl1*(x3-x2)-xl2*(y3-y2); y=y3+xl2*(x3-x2)+xl1*(y3-y2); // Test for convergence if (fabs(x-x0)+fabs(y-yy0)=n) return; // Continue k=k+1; x0=x; yy0=y; x1=x2; y1=y2; x2=x3; y2=y3; x3=x; y3=y; goto e100; } // ZMueller() void main() { printf("\n Write the initial guesses and their bounds:\n\n"); printf(" X0 = "); scanf("%lf",&x0); printf(" Bond on X0 = "); scanf("%lf",&b1); printf("\n"); printf(" yy0 = "); scanf("%lf",&yy0); printf(" Bond on yy0 = "); scanf("%lf",&b2); printf("\n Convergence criterion: "); scanf("%lf",&e); printf("\n Maximum number of iterations: "); scanf("%d",&n); ZMueller(); // Call complex domain Mueller subroutine printf("\n\n The calculated zero is (X,Y) =\n"); printf(" %f %f\n\n",x,y); Z(x,y,&u,&v); printf(" The associated Z value is (U,V) =\n"); printf(" %f %f\n\n",u,v); printf(" The number of steps was: %d\n\n",k); } // End of file zmueller.cpp