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 ******