/****************************************************** * Program to demonstrate the zero searching algorithm * * --------------------------------------------------- * * 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) * * --------------------------------------------------- * * SAMPLE RUN: * * What is the initial guess for x and y (x0,y0): * * * * X0 = 0 * * Y0 = 0 * * * * What is the radius of the first search circle: 5 * * * * By what fraction this circle will be reduced on * * each iteration: .5 * * * * How many evaluation points per quadrant: 4 * * * * Maximum number of iterations: 10 * * * * * * The approximate solution is: * * Z = 0.000576 + 0.999857 I * * * * Number of iterations: 10 * * * ******************************************************/ #include // for printf(), scanf() #include // for strcpy() #include // for cos(), sin() // global variables double e,w,x0,yy0,xx,yy; int k,m,n; char z[4]; //*************************************************** // Complex function(x,y) subroutine void Eval(double x, double y, double *u, double *v) { *u=x*x-y*y+1.0; *v=2.0*x*y; } //*************************************************** // Local functions used by Zcircle() void Set_to_zero(double *x,double *y) { *x=0.0; *y=0.0; } void Auxilliary(int m, int *i,int *j) { *i=*i+1; *j=*j-1; if (*i>m) *i=*i-m; if (*j<1) *j=*j+m; } void Interpolation(int m, int i,int j,double *X,double *Y, double *U,double *V,double *xx,double *yy) { j=i-1; if (j<1) j=j+m; //Regula Falsi interpolation for the zero *xx=(X[i]*U[j]+X[j]*U[i])/(U[i]+U[j]); *yy=(Y[i]*V[j]+Y[j]*V[i])/(V[i]+V[j]); } // End local functions /************************************************ * Complex root search subroutine * * --------------------------------------------- * * This routine searches for the complex roots * * of an analytical function by encircling the * * zero and estimating where it is. The circle * * is subsequently tightened by a factor e, and * * a new estimate made. * * The inputs to the subroutine are; * * w initial radius of search circle * * x0,y0 the initial guess * * e factor by which circle is reduced * * n maximum number of iterations * * m evaluation points per quadrant * * The routine returns Z=X+IY (X,Y), and the * * number of iterations, k. * ************************************************/ void Zcircle() { //Labels: 50,100,200,250,300,310,320,330,340,350,400,450 // 460,470,500,510,520,521 int i,j,m3,m4; double X[80],Y[80],U[80],V[80]; double a,b1,b2,PI,x1,x2,xm1,xm2; double x3,y3,x4,y1,y2,y4,uu,vv; m=m*4; k=1; PI = 4*atan(1); a=2*PI/m; e50:for (i=1; im) i=i-m; j=i; uu=U[i]; // Flip directions alternately e250: Auxilliary(m,&i,&j); if (uu*U[i]<0) goto e300; if (uu*U[j]<0) goto e310; // Guard against infinite loop if (i==xm1+(xm2/2)) Set_to_zero(&xx,&yy); if (j==xm1+(xm2/2)) Set_to_zero(&xx,&yy); goto e250; // Transition found e300: m3=i; goto e320; // Transition found e310: if (j==m) j=0; m3=j+1; // xm1 and m3 have been determined, now for xm2 and m4 // now for the vv transitions e320: i=(int) xm1+m/4; if (i>m) i=i-m; j=i; vv=V[i]; e330: Auxilliary(m,&i,&j); if (vv*V[i]<0) goto e340; if (vv*V[j]<0) goto e350; // Guard against infinite loop if (i==xm1+m/4) Set_to_zero(&xx,&yy); if (j==xm1+m/4) Set_to_zero(&xx,&yy); goto e330; // xm2 has been found e340: xm2=i; goto e400; e350: if (j==m) j=0; xm2=j+1; // xm2 has been found, now for m4 e400: i=(int) xm2+m/2; if (i>m) i=i-m; j=i; vv=V[i]; e450: Auxilliary(m,&i,&j); if (uu*V[i]<0) goto e460; if (vv*V[j]<0) goto e470; // Guard against infinite loop if (i==xm2+m/2) Set_to_zero(&xx,&yy); if (j==xm2+m/2) Set_to_zero(&xx,&yy); goto e450; e460: m4=i; goto e500; e470: if (j==m) j=0; m4=j+1; // All the intersections have been determine // Interpolate to find the four (x,y) coordinates e500: i=(int) xm1; Interpolation(m,i,j,X,Y,U,V,&xx,&yy); x1=xx ; y1=yy ; i=(int) xm2; Interpolation(m,i,j,X,Y,U,V,&xx,&yy); x2=xx ; y2=yy ; i=m3; Interpolation(m,i,j,X,Y,U,V,&xx,&yy); x3=xx ; y3=yy ; i=m4; Interpolation(m,i,j,X,Y,U,V,&xx,&yy); x4=xx ; y4=yy; // Calculate the intersection of the lines // Guard against a divide by zero if (x1!=x3) goto e510; xx=x1; yy=(y1+y3)/2.0; goto e520; e510: xm1=(y3-y1)/(x3-x1); if (x2!=x4) goto e520; xm2=1e8; goto e521; e520: xm2=(y2-y4)/(x2-x4); e521: b1=y1-xm1*x1; b2=y2-xm2*x2; xx=-(b1-b2)/(xm1-xm2); yy=(xm1*b2+xm2*b1)/(xm1+xm2); // is another iteration in order ? if (k==n) return; x0=xx ; yy0=yy ; k=k+1 ; w=w*e; goto e50; } // Zcircle void main() { printf("\n What is the initial guess for x and y (x0,yy0):\n\n"); printf(" X0 = "); scanf("%lf",&x0); printf(" Y0 = "); scanf("%lf",&yy0); printf("\n What is the radius of the first search circle: "); scanf("%lf",&w); printf("\n By what fraction this circle will be reduced on each iteration: "); scanf("%lf",&e); printf("\n How many evaluation points per quadrant: "); scanf("%d",&m); printf("\n Maximum number of iterations: "); scanf("%d",&n); printf("\n"); Zcircle(); // Call ZCircle function printf("\n"); if (yy>=0) strcpy(z," + "); if (yy< 0) strcpy(z," - "); printf("\n The approximate solution is:\n"); printf(" Z = %f%s%f I\n\n",xx,z,yy); printf(" Number of iterations: %d\n\n",k); } // End of file Zcircle.cpp