/****************************************************** * Program to demonstrate the Newton root subroutine * * in the complex domain * * --------------------------------------------------- * * 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: * * * * (Example: find a complex root of z^2 + 1 = 0) * * * * What is the initial guess: * * * * X0 = .2 * * Y0 = 1.2 * * * * Convergence criterion: 1e-8 * * Maximum number of iterations: 10 * * * * The root estimate is: * * * * X0 = -0.000000 * * Y0 = 1.000000 * * * * The number of iterations performed was: 5 * * * ******************************************************/ #include #include #define TINY 1e-12 //small number int k,n; double a,e,x0,yy0,u,v,u1,v1,u2,v2,x,y; //********************************************************* // Functions subroutine void Eval(double x, double y, double *u, double *v, double *u1, double *v1, double *u2, double *v2) { *u=x*x-y*y+1.0; *v=2.0*x*y; *u1=2.0*x; *u2=-2.0*y; *v1=2.0*y; *v2=2.0*x; } //********************************************************* /************************************************ * Complex root seeking using Newton"s method * * --------------------------------------------- * * This routine uses the complex domain form of * * Newton"s method for iteratively searching * * for roots. The complex function and its first * * partial derivatives must be available in the * * form; F(Z) = U(X,Y) + I V(X,Y). The required * * derivatives are DU/DX and DU/DY. * * --------------------------------------------- * * INPUTS; initial guess, x0 and y0, convergence * * criteria, e, maximum number of iterations, n. * * OUTPUTS; approximation to the root, x and y, * * number of performed iterations, k. * ************************************************/ void ZNewton() { // Label: e100 k=0; e100: k++; // Get u,v and the derivatives u1,v1,u2,v2 Eval(x0,yy0,&u,&v,&u1,&v1,&u2,&v2); a=u1*u1+u2*u2; // Guard against a=0 if (a < TINY) { printf("\n ZERO DIVIDE ERROR - u1*u1+u2*u2 must be <> 0.\n\n"); return; } x=x0+(v*u2-u*u1)/a; y=yy0-(v*u1+u*u2)/a; // Check for convergence in euclidean space if ((x0-x)*(x0-x)+(yy0-y)*(yy0-y)<=e*e) return; if (k>=n) return; x0=x ; yy0=y; goto e100; } // ZNewton() void main() { printf("\n What is the initial guess:\n\n"); printf(" X0 = "); scanf("%lf",&x0); printf(" Y0 = "); scanf("%lf",&yy0); printf("\n Convergence criterion: "); scanf("%lf",&e); printf("\n Maximum number of iterations: "); scanf("%d",&n); ZNewton(); // Call ZNewton subroutine if (a!=0) { printf("\n\n The root estimate is:\n\n"); printf(" X0 = %f\n",x); printf(" Y0 = %f\n",y); printf("\n The number of iterations performed was: %d\n\n",k); } } // End of file Znewton.cpp