!**************************************************** !* 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.259967402504241 -0.399661390631072 * !* -0.111845742991334 0.668006959856486 * !* -1.07915614227153 4.90406822387803 * !* * !* ------------------------------------------------ * !* Ref.: From Numath Library By Tuan Dang Trong in * !* Fortran 77 [BIBLI 18]. * !* * !* F90 Release By J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !**************************************************** PROGRAM TEST_NEWTON parameter(N=3) real*8 A(N+1), B(N+1), RR(N+1), RI(N+1) real*8 WORK(4*N+5) integer IERR A(1)= 2.5d0; B(1)=1.d0 A(2)= 7.5d0; B(2)=-12.d0 A(3)=-3.75d0; B(3)=0.4d0 A(4)= 4.25d0; B(4)=-1.d0 ! Example #2 (real roots are -3, 1, 2) ! A(1)= 1.d0; B(1)=0.d0 ! A(2)= 0.d0; B(2)=0.d0 ! A(3)=-7.d0; B(3)=0.d0 ! A(4)= 6.d0; B(4)=0.d0 call NEWTON (N,A,B,RR,RI,IERR,WORK) print *,' ' print *,' Error code =', IERR print *,' ' print *,' Complex roots are:' do I=1, N write(*,*) RR(I), RI(I) end do print *,' ' stop END SUBROUTINE NEWTON (N,A,B,X,Y,IER,W) IMPLICIT REAL *8 (A-H,O-Z) !--------------------------------------------------------------------- ! 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 ZONE: ! W TABLE OF DIMENSION 4*N+5 ! REFERENCE: ! E.DURAND -SOLUTIONS NUMERIQUES DES EQUATIONS ALGEBRIQUES, ! TOME 1, MASSON & CIE, 1960, PP.260-266 !---------------------------------------------------------------------- DIMENSION A(*),B(*),X(*),Y(*),W(*) N1 = N+2 N2 = N+N+3 N3 = N+N+N+4 CALL CNWTON (A,B,X,Y,N,IER,W(1),W(N1),W(N2),W(N3)) RETURN END ! SUBROUTINE CNWTON (A,B,X,Y,N,IER,G,H,T,F) IMPLICIT REAL *8 (A-H,O-Z) DIMENSION A(*),B(*),X(*),Y(*),G(*),H(*),T(*),F(*) DATA IT1 /25/,IT2 /100/,TOL /1.D-2/,EPS /0.D0/ IER = 0 IF (EPS.NE.0.D0) GO TO 7 EPS = 1.D0 5 EPS = 0.5D0*EPS EPS1 = EPS+1.D0 IF (EPS1.GT.1.D0) GO TO 5 7 L = N X0 = 1.D0 Y0 = 1.D0 G(1) = A(1) T(1) = A(1) H(1) = B(1) F(1) = B(1) 10 IF (L.EQ.1) GO TO 55 M = L+1 I = N-L+1 K1 = 0 K2 = 0 X1 = X0 Y1 = Y0 15 DO 20 K = 2,M G(K) = A(K)+X1*G(K-1)-Y1*H(K-1) 20 H(K) = B(K)+Y1*G(K-1)+X1*H(K-1) DO 25 K = 2,L T(K) = G(K)+X1*T(K-1)-Y1*F(K-1) 25 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 = (DABS(DX)+DABS(DY))/(DABS(X2)+DABS(Y2)) K2 = K2+1 X1 = X2 Y1 = Y2 IF (R1.LT.TOL) GO TO 30 IF (K2.LE.IT2) GO TO 15 IER = I GO TO 60 30 K1 = K1+1 IF (K1.LE.IT1) GO TO 15 X(I) = X2 Y(I) = Y2 IF (DABS(X2).LE.EPS) X2 = EPS R2 = DABS(Y2/X2) IF (R2.LT.0.1D0) GO TO 40 IF (R2*EPS.GT.1.D0) X(I) = 0.D0 X0 = X2 Y0 = Y2 GO TO 45 40 IF (R2.LE.EPS) Y(I) = 0.D0 X0 = X2+Y2 Y0 = X2-Y2 45 DO 50 K = 1,L A(K) = G(K) 50 B(K) = H(K) L = L-1 GO TO 10 55 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 (DABS(X(N)).LE.EPS) X(N) = EPS R2 = DABS(Y(N)/X(N)) IF (R2*EPS.GT.1.D0) X(N) = 0.D0 IF (R2.LE.EPS) Y(N) = 0.D0 60 CONTINUE DO L=1,N P=X(L) Q=Y(L) R=P*P+Q*Q IF(L.EQ.1) GO TO 70 DO 65 I=L,2,-1 CABS=X(I-1)**2+Y(I-1)**2 IF(R.GE.CABS) GO TO 75 X(I)=X(I-1) Y(I)=Y(I-1) 65 CONTINUE 70 I=1 75 X(I)=P Y(I)=Q ENDDO END !end of file tnewton.f90