home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
BURKS 2
/
BURKS_AUG97.ISO
/
BURKS
/
LANGUAGE
/
FORTRAN
/
F77TO90
/
code
/
lab2_e.f
< prev
next >
Wrap
Text File
|
1996-01-09
|
6KB
|
116 lines
******************************************************************************
** 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 ******