/**************************************************** * Program to demonstrate the complex domain * * Allroot 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 the two complex roots of zē+1 = 0 * * * * SAMPLE RUN: * * * * Input the initial guesses and their bounds: * * * * X0 = 4 * * Bond on X0 = 1 * * Y0 = 4 * * Bond on Y0 = 1 * * * * Convergence criterion: 1e-6 * * Maximum number of iterations: 30 * * How many roots are to be sought: 2 * * Is the function defined in 1000 (1) * * or is it a series (2): 1 * * * * The estimated roots are: * * * * X = -0.000018 * * Y = 1.000032 * * * * X = 4.020035 * * Y = 4.028881 * * * * (Note: only one root is found) * ****************************************************/ #include #include #define PI 4*atan(1) double A[10], X[10], Y[10]; double b1,b2,e,u,v,u1,v1,u2,v2,u3,v3,x0,y00,xx,yy,x4,y4,z1,z2; int i,jj1,k,m,n,n2,n3; // Rectangular to polar conversion void RectPol(double x, double y, double *u, double *v) { *u=sqrt(x*x+y*y); //Guard against ambiguous vector if (y==0.0) y=1e-16; //Guard against divide by zero if (x==0.0) x=1e-16; *v=atan(y/x); //Check quadrant and adjust if (x<0.0) *v=*v+PI; if (*v<0.0) *v=*v+2.0*PI; } //Polar to rectangular conversion void PolRect(double u, double v, double *x, double *y) { *x=u*cos(v); *y=u*sin(v); } //Polar power void PolPower(int n, double u, double v, double *u1, double *v1) { int i; *u1=u; for (i=2; ia2) goto e200; a1=e1-b1; a2=e2-b2; goto e300; e200: a1=e1+b1; a2=e2+b2; e300: aa=a1*a1+a2*a2; xl1=a1*d1*u3-a1*d2*v3+a2*u3*d2+a2*v3*d1; //Guard against divide by zero if (aa==0.0) aa=1e-7; xl1=-2.0*xl1/aa; xl2=-d1*u3*a2+d2*v3*a2+a1*u3*d2+a1*v3*d1; xl2=-2.0*xl2/aa; //Calculate new estimate xx=x3+xl1*(x3-x2)-xl2*(y3-y2); yy=y3+xl2*(x3-x2)+xl1*(y3-y2); //Test for convergence if (fabs(xx-x0)+fabs(yy-y00)=n) return; //Continue k=k+1; x0=xx; y00=yy; x1=x2; y1=y2; x2=x3; y2=y3; x3=xx; y3=yy; goto e100; } //Zmueller //***** Coefficients subroutine ***** void Coefficients() { m=5; A[0]=0; A[1]=24; A[2]=-50; A[3]=35; A[4]=-10; A[5]=1; } /*************************************************** * General root determination subroutine * * ------------------------------------------------ * * This routine attempts to calculate the several * * roots of a given series or function by repea- * * tedly using the Zmueller subroutine and removing * * the roots already found by division. * * ------------------------------------------------ * * INPUTS: X0, Y0 - The initial guess * * b1, b2 - The bounds on this guess * * e - The convergence criteria * * n - The maximum number of iterations * * for each root * * n2 - the number of roots to seek * * n3 - A flag (1) for a function F(Z) * * (2) for a polynomial. * * ------------------------------------------------ * * OUTPUTS: The n2 roots found X(i), Y(i) * * The last number of iterations, k. * ***************************************************/ void AllRoot() { //Labels: e100, e200; k=0; if (n3==1) goto e100; if (n3!=2) return; //error in n3 e100: jj1=0; //Save the initial guess x4=x0; y4=y00; //If n3=2 get the series coefficients if (n3==2) Coefficients(); //Test for completion e200: if (jj1==n2) return; //Call Zmueller subroutine ZMueller(); jj1++; X[jj1]=xx; Y[jj1]=yy; x0=x4; y00=y4; //Try another pass goto e200; } 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(" Y0 = "); scanf("%lf",&y00); printf(" Bond on Y0 = "); scanf("%lf",&b2); printf("\n Convergence criterion: "); scanf("%lf",&e); printf(" Maximum number of iterations: "); scanf("%d",&n); printf(" How many roots are to be sought: "); scanf("%d",&n2); printf("\n Is the function defined in Eval (1)"); printf(" or is it a series (2) "); scanf("%d",&n3); AllRoot(); //Call complex domain Allroot subroutine printf("\n\n The estimated roots are:"); for (i=1; i