home *** CD-ROM | disk | FTP | other *** search
/ BURKS 2 / BURKS_AUG97.ISO / BURKS / LANGUAGE / FORTRAN / F77TO90 / code / lab2_e.f < prev    next >
Text File  |  1996-01-09  |  6KB  |  116 lines
  1. ******************************************************************************
  2. **      This program calculates all roots (real and/or complex) to an       **
  3. **      N:th-degree polynomial with complex coefficients. (N <= 10)         **
  4. **                                                                          **
  5. **                 n        n-1                                             **
  6. **       P(z) = a z  + a   z    + ... + a z + a                             **
  7. **               n      n-1              1     0                            **
  8. **                                                                          **
  9. **      The execution terminates if                                         **
  10. **                                                                          **
  11. **                 1) Abs (Z1-Z0) < EPS       ==>   Root found = Z1         **
  12. **                 2) ITER > MAX              ==>   Slow convergence        **
  13. **                                                                          **
  14. **      The program sets EPS to 1.0E-7 and MAX to 30                        **
  15. **                                                                          **
  16. **      The NEWTON-RAPHSON method is used:   z1 = z0 - P(z0)/P'(z0)         **
  17. **      The values of P(z0) and P'(z0) are calculated with HORNER'S SCHEME. **
  18. **                                                                          **
  19. **      The array A(0:10) contains the complex coefficients of the          ** 
  20. **      polynomial P(z).                                                    **
  21. **      The array B(1:10) contains the complex coefficients of the          **
  22. **      polynomial Q(z),                                                    **
  23. **                  where P(Z) = (z-z0)*Q(z) + P(z0)                        **
  24. **      The coefficients to Q(z) are calculated with HORNER'S SCHEME.       **
  25. **                                                                          **
  26. **      When the first root has been obtained with the NEWTON-RAPHSON         **
  27. **      method, it is divided away (deflation).                             **
  28. **      The quotient polynomial is Q(z).                                    **
  29. **      The process is repeated, now using the coefficients of Q(z) as      **
  30. **      input data.                                                         **
  31. **      As a starting approximation the value STARTV = 1+i is used          **
  32. **      in all cases.                                                       **
  33. **      Z0 is the previous approximation to the root.                       **
  34. **      Z1 is the latest approximation to the root.                         **
  35. **      F0 = P(Z0)                                                          **
  36. **      FPRIM0 = P'(Z0)                                                     **
  37. ******************************************************************************
  38.       COMPLEX      A(0:10), B(1:10), Z0, Z1, STARTV
  39.       INTEGER      N, I, ITER, MAX
  40.       REAL         EPS
  41.       DATA  EPS/1E-7/, MAX /30/, STARTV /(1,1)/
  42. ******************************************************************************
  43. 20    WRITE(6,*)  'Give the degree of the polynomial'
  44.       READ (5,*) N
  45.       IF (N .GT. 10) THEN
  46.          WRITE(6,*) 'The polynomial degree is maximum 10'
  47.       GOTO 20
  48.       WRITE (6,*) 'Give the coefficients of the polynomial,',
  49.      , ' as complex constants'
  50.       WRITE (6,*) 'Highest degree coefficient first'
  51.       DO 30 I = N, 0, -1
  52.           WRITE (6,*) 'A(' , I, ') = '
  53.           READ  (5,*)  A(I)
  54. 30    CONTINUE
  55.       WRITE (5,*) '    The roots are','        Number of iterations'
  56. ******************************************************************************
  57. 40    IF (N GT 0) THEN
  58. C     ******     Find the next root                                     ******
  59.           Z0 = (0,0)
  60.           ITER = 0
  61.           Z1 = STARTV
  62. 50        IF (ABS(Z1-Z0) .GE. EPS) THEN 
  63. C         ++++++ Continue the iteration until two estimates             ++++++
  64. C         ++++++ are sufficiently close.                                ++++++
  65.               ITER = ITER + 1
  66.               IF (ITER .GT. MAX) THEN
  67. C             ------ Too many iterations  ==> TERMINATE                 ------
  68.                      WRITE (6,*) 'Too many iterations.'
  69.                      WRITE (6,*) 'The latest approximation to the root is ',Z1
  70.               GOTO 200
  71.               ENDIF
  72.               Z0 = Z1
  73.               HORNER (N, A, B, Z0, F0, FPRIM0)
  74. C             ++++++   NEWTON-RAPHSON'S METHOD                          ++++++
  75.               Z1 = Z0 - F0/FPRIM0
  76.           GOTO 50
  77.           ENDIF
  78. C         ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  79. 100       WRITE (6,*) Z1, '         ',Iter
  80. C         ****** The root is found. Divide it away and seek the next one *****
  81.           N = N - 1
  82.           FOR I = 0 TO N DO
  83.               A(I) = B(I+1)
  84.       GOTO 40
  85.       ENDIF
  86. 200   END
  87.  
  88.       SUBROUTINE  HORNER(N, A, B, Z, F, FPRIM)
  89. ******    For the meaning of the parameters - see the main program      ******
  90. ******    BI and CI are temporary variable.                             ******
  91.  
  92.       INTEGER N, I
  93.       COMPLEX A(1:10), B(0:10), Z, F, FPRIM, BI, CI
  94.  
  95.       BI = A(N)
  96.       B(N) = BI
  97.       CI = BI
  98.  
  99.       DO 60  I = N-1, 1, -1
  100.              BI = A(I) + Z*BI
  101.              CI = BI + Z*CI
  102.              B(I) = BI
  103. C            ++++++ BI = B(I) for the calculation of P(Z)               ++++++
  104. C            ++++++ CI  for the calculation of P'(Z)                    ++++++
  105. 60    CONTINUE
  106.  
  107.       F = A(0) + Z*BI
  108.       FPRIM = CI
  109.  
  110.       RETURN
  111.       END
  112. ******   END HORNER'S SCHEME                                            ******
  113.  
  114. ******   This program is composed by Ulla Ouchterlony 1984              ******
  115. ******   The program is translated by Bo Einarsson 1995                 ******
  116.