PROGRAM LAB2 !***************************************************************************** !* This program contains some intentional errors. ** !* ** !* 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 ! ****** Find the next root ****** Z0 = (0,0) ITER = 0 Z1 = STARTV 50 IF (ABS(Z1-Z0) .GE. EPS) THEN ! ++++++ Continue the iteration until two estimates ++++++ ! ++++++ are sufficiently close. ++++++ ITER = ITER + 1 IF (ITER .GT. MAX) THEN ! ------ 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) ! ++++++ NEWTON-RAPHSON'S METHOD ++++++ Z1 = Z0 - F0/FPRIM0 GOTO 50 ENDIF ! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 100 WRITE (6,*) Z1, ' ',Iter ! ****** 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 PROGRAM LAB2 SUBROUTINE HORNER(N, A, B, Z, F, FPRIM) !***** For the meaning of the parameters - see the main program ****** !***** BI and CI are temporary variables. ****** 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 ! ++++++ BI = B(I) for the calculation of P(Z) ++++++ ! ++++++ CI for the calculation of P'(Z) ++++++ 60 CONTINUE F = A(0) + Z*BI FPRIM = CI RETURN END SUBROUTINE HORNER !***** END HORNER'S SCHEME ****** !***** This program was composed by Ulla Ouchterlony in 1984 ****** !***** The program was translated by Bo Einarsson in 1995 ******