/*************************************************** * Find all roots of a complex polynomial using * * Newton's iterative formulation. * * ------------------------------------------------ * * SAMPLE RUN: * * (Find roots of complex polynomial: * * (2.5 + I) X3 + (7.5 - 12 I) X2 + (-3.75 + 0.4 I * * ) X + (4.25 - I) ) * * * * Error code = 0 * * * * Complex roots are: * * 0.25996740 -0.39966139 * * -0.11184574 0.66800696 * * -1.07915614 4.90406822 * * * * ------------------------------------------------ * * Ref.: From Numath Library By Tuan Dang Trong in * * Fortran 77 [BIBLI 18]. * * * * C++ Release By J-P Moreau, Paris. * * (www.jpmoreau.fr) * ***************************************************/ #include #include #define NMAX 50 int I,IERR, N; //IERR: error code, N: order of polynomial double A[NMAX],B[NMAX],RR[NMAX],RI[NMAX]; /* A: real parts of complex coefficients B: imaginary parts of complex coefficients (stored in decreasing powers) RR: real parts of roots RI: imaginary parts of roots */ double W1[NMAX],W2[NMAX],W3[NMAX],W4[NMAX]; //working zones void CNEWTON (int N, double *A, double *B, double *X, double *Y, int *IER, double *G, double *H, double *T, double *F) { /*-------------------------------------------------------------------- ! CALCULATE ALL THE ROOTS OF A COMPLEX POLYNOMIAL USING ! NEWTON'S ITERATIVE FORMULATION ! ! INPUTS: ! N ORDER OF POLYNOMIAL ! A,B TABLES OF DIMENSION N+1, RESPECTIVELY REAL, ! IMAGINARY PART OF POLYNOMIAL COMPLEX COEFFICIENTS ! C = A+I*B, STORED IN DECREASING POWERS. ! OUTPUTS: ! X,Y TABLES OF DIMENSION N, RESPECTIVELY REAL, ! IMAGINARY PART OF COMPLEX ROOTS FOUND, Z = X+I*Y ! IER ERROR CODE = 0, CONVERGENCE OK ! = I, NO CONVERGENCE FOR I-TH ROOT ! WORKING ZONES: ! G,H,T,F TABLES OF DIMENSION N !---------------------------------------------------------------------- ! REFERENCE: ! E.DURAND -SOLUTIONS NUMERIQUES DES EQUATIONS ALGEBRIQUES, ! TOME 1, MASSON & CIE, 1960, PP.260-266 !---------------------------------------------------------------------*/ // Labels: e5,e7,e10,e15,e30,e40,e45,e55,e60,e70,e75 int I,IT1,IT2,K,K1,K2,L,M; double EPS,EPS1,P,Q,R,TOL,X0,X1,Y0,Y1; double CABS,DX,DY,GD,GX,GY,R1,R2,R3,X2,Y2; IT1=25; IT2=100; TOL=1e-2; EPS=0.0; IER = 0; if (EPS != 0.0) goto e7; EPS = 1.0; e5: EPS = 0.50*EPS; EPS1 = EPS+1.0; if (EPS1 > 1.0) goto e5; e7: L = N; X0 = 1.0; Y0 = 1.0; G[1] = A[1]; T[1] = A[1]; H[1] = B[1]; F[1] = B[1]; e10: if (L == 1) goto e55; M = L+1; I = N-L+1; K1 = 0; K2 = 0; X1 = X0; Y1 = Y0; e15: for (K=2; K<=M; K++) { G[K] = A[K]+X1*G[K-1]-Y1*H[K-1]; H[K] = B[K]+Y1*G[K-1]+X1*H[K-1]; } for (K=2; K<=L; K++) { T[K] = G[K]+X1*T[K-1]-Y1*F[K-1]; F[K] = H[K]+Y1*T[K-1]+X1*F[K-1]; } GD = T[L]*T[L]+F[L]*F[L]; GX = G[M]*T[L]+H[M]*F[L]; GY = H[M]*T[L]-G[M]*F[L]; DX = -GX/GD; DY = -GY/GD; X2 = X1+DX; Y2 = Y1+DY; R1 = (fabs(DX)+fabs(DY))/(fabs(X2)+fabs(Y2)); K2 = K2+1; X1 = X2; Y1 = Y2; if (R1 < TOL) goto e30; if (K2 <= IT2) goto e15; *IER = I; goto e60; e30: K1++; if (K1 <= IT1) goto e15; X[I] = X2; Y[I] = Y2; if (fabs(X2) <= EPS) X2 = EPS; R2 = fabs(Y2/X2); if (R2 < 0.1) goto e40; if (R2*EPS > 1.0) X[I] = 0.0; X0 = X2; Y0 = Y2; goto e45; e40: if (R2 <= EPS) Y[I] = 0.0; X0 = X2+Y2; Y0 = X2-Y2; e45: for (K=1; K<=L; K++) { A[K] = G[K]; B[K] = H[K]; } L--; goto e10; e55: R3 = G[1]*G[1]+H[1]*H[1]; X[N] = -(G[1]*G[2]+H[1]*H[2])/R3; Y[N] = (H[1]*G[2]-G[1]*H[2])/R3; if (fabs(X[N]) <= EPS) X[N] = EPS; R2 = fabs(Y[N]/X[N]); if (R2*EPS > 1.0) X[N] = 0.0; if (R2 <= EPS) Y[N] = 0.0; e60: for (L=1; L<=N; L++) { P=X[L]; Q=Y[L]; R=P*P+Q*Q; if (L == 1) goto e70; for (I=L; I>1; I--) { CABS=X[I-1]*X[I-1]+Y[I-1]*Y[I-1]; if (R >= CABS) goto e75; X[I]=X[I-1]; Y[I]=Y[I-1]; } e70: I=1; e75: X[I]=P; Y[I]=Q; } } //CNEWTON void main() { N = 3; //order of polynomial //define vectors A, B of coefficients A[1]= 2.5; B[1]=1.0; A[2]= 7.5; B[2]=-12.0; A[3]=-3.75; B[3]=0.4; A[4]= 4.25; B[4]=-1.0; /* Example #2 (real roots are -3, 1, 2) A[1]= 1.; B[1]=0.0; A[2]= 0.; B[2]=0.0; A[3]=-7.; B[3]=0.0; A[4]= 6.; B[4]=0.0; */ //call main function CNEWTON(N,A,B,RR,RI,&IERR,W1,W2,W3,W4); //print results printf("\n Error code = %d\n", IERR); printf("\n Complex roots are:\n"); for (I=1; I<=N; I++) printf(" %10.8f %10.8f\n", RR[I], RI[I]); printf("\n"); } // end of file tnewton.cpp