****************************************************************************** ** This program calculates all roots (real and/or complex) to an ** ** N:th-degree polynomial with complex coefficients. (N <= 10) ** ** ** ** n n-1 ** ** P(z) = a z + a z + ... + a z + a ** ** n n-1 1 0 ** ** ** ** The execution terminates if ** ** ** ** 1) Abs (Z1-Z0) < EPS ==> Root found = Z1 ** ** 2) ITER > MAX ==> Slow convergence ** ** ** ** The program sets EPS to 1.0E-7 and MAX to 30 ** ** ** ** The NEWTON-RAPHSON method is used: z1 = z0 - P(z0)/P'(z0) ** ** The values of P(z0) and P'(z0) are calculated with HORNER'S SCHEME. ** ** ** ** The array A(0:10) contains the complex coefficients of the ** ** polynomial P(z). ** ** The array B(1:10) contains the complex coefficients of the ** ** polynomial Q(z), ** ** where P(Z) = (z-z0)*Q(z) + P(z0) ** ** The coefficients to Q(z) are calculated with HORNER'S SCHEME. ** ** ** ** When the first root has been obtained with the NEWTON-RAPHSON ** ** method, it is divided away (deflation). ** ** The quotient polynomial is Q(z). ** ** The process is repeated, now using the coefficients of Q(z) as ** ** input data. ** ** As a starting approximation the value STARTV = 1+i is used ** ** in all cases. ** ** Z0 is the previous approximation to the root. ** ** Z1 is the latest approximation to the root. ** ** F0 = P(Z0) ** ** FPRIM0 = P'(Z0) ** ****************************************************************************** COMPLEX A(0:10), B(1:10), Z0, Z1, STARTV INTEGER N, I, ITER, MAX REAL EPS DATA EPS/1E-7/, MAX /30/, STARTV /(1,1)/ ****************************************************************************** 20 WRITE(6,*) 'Give the degree of the polynomial' READ (5,*) N IF (N .GT. 10) THEN WRITE(6,*) 'The polynomial degree is maximum 10' GOTO 20 WRITE (6,*) 'Give the coefficients of the polynomial,', , ' as complex constants' WRITE (6,*) 'Highest degree coefficient first' DO 30 I = N, 0, -1 WRITE (6,*) 'A(' , I, ') = ' READ (5,*) A(I) 30 CONTINUE WRITE (5,*) ' The roots are',' Number of iterations' ****************************************************************************** 40 IF (N GT 0) THEN C ****** Find the next root ****** Z0 = (0,0) ITER = 0 Z1 = STARTV 50 IF (ABS(Z1-Z0) .GE. EPS) THEN C ++++++ Continue the iteration until two estimates ++++++ C ++++++ are sufficiently close. ++++++ ITER = ITER + 1 IF (ITER .GT. MAX) THEN C ------ Too many iterations ==> TERMINATE ------ WRITE (6,*) 'Too many iterations.' WRITE (6,*) 'The latest approximation to the root is ',Z1 GOTO 200 ENDIF Z0 = Z1 HORNER (N, A, B, Z0, F0, FPRIM0) C ++++++ NEWTON-RAPHSON'S METHOD ++++++ Z1 = Z0 - F0/FPRIM0 GOTO 50 ENDIF C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 100 WRITE (6,*) Z1, ' ',Iter C ****** The root is found. Divide it away and seek the next one ***** N = N - 1 FOR I = 0 TO N DO A(I) = B(I+1) GOTO 40 ENDIF 200 END SUBROUTINE HORNER(N, A, B, Z, F, FPRIM) ****** For the meaning of the parameters - see the main program ****** ****** BI and CI are temporary variable. ****** INTEGER N, I COMPLEX A(1:10), B(0:10), Z, F, FPRIM, BI, CI BI = A(N) B(N) = BI CI = BI DO 60 I = N-1, 1, -1 BI = A(I) + Z*BI CI = BI + Z*CI B(I) = BI C ++++++ BI = B(I) for the calculation of P(Z) ++++++ C ++++++ CI for the calculation of P'(Z) ++++++ 60 CONTINUE F = A(0) + Z*BI FPRIM = CI RETURN END ****** END HORNER'S SCHEME ****** ****** This program is composed by Ulla Ouchterlony 1984 ****** ****** The program is translated by Bo Einarsson 1995 ******