home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / villad / jcobi.f < prev    next >
Text File  |  1996-09-28  |  6KB  |  241 lines
  1. ****************************************************************
  2. *
  3. *     The following routines (JCOBI, DIF, DFOPR, INTRP, AND RADAU)
  4. *     are the same as found in Villadsen, J. and M.L. Michelsen,
  5. *     Solution of Differential Equation Models by Polynomial
  6. *     Approximation, Prentice-Hall (1978) pages 418-420.
  7. *
  8. *     Cosmetic changes (elimination of arithmetic IF statements, most
  9. *     GO TO statements, and indentation of program blocks) made by:
  10. *
  11. *     John W. Eaton
  12. *     Department of Chemical Engineering
  13. *     The University of Texas at Austin
  14. *     Austin, Texas 78712
  15. *
  16. *     June 6, 1987
  17. *
  18. *     Some error checking additions also made on June 7, 1987
  19. *
  20. *     Further cosmetic changes made August 20, 1987
  21. *
  22. ************************************************************************
  23. *
  24.       SUBROUTINE JCOBI
  25.      +  (
  26.      +  ND, N, N0, N1, ALPHA, BETA, DIF1, DIF2, DIF3, ROOT
  27.      +  )
  28. C
  29.       INTEGER
  30.      +
  31.      +  ND, N, N0, N1
  32. C 
  33.       DOUBLE PRECISION
  34.      +
  35.      +  ALPHA, BETA, DIF1(ND), DIF2(ND), DIF3(ND), ROOT(ND)
  36. C
  37. C***********************************************************************
  38. C
  39. C     VILLADSEN AND MICHELSEN, PAGES 131-132, 418
  40. C
  41. C     THIS SUBROUTINE COMPUTES THE ZEROS OF THE JACOBI POLYNOMIAL
  42. C
  43. C        (ALPHA,BETA)
  44. C       P  (X)
  45. C        N
  46. C
  47. C     USE DIF (GIVEN BELOW) TO COMPUTE THE DERIVATIVES OF THE NODE
  48. C     POLYNOMIAL
  49. C
  50. C                     N0     (ALPHA,BETA)           N1
  51. C       P  (X)  =  (X)   *  P (X)         *  (1 - X)
  52. C        NT                   N
  53. C
  54. C     AT THE INTERPOLATION POINTS.
  55. C
  56. C     INPUT PARAMETERS:
  57. C
  58. C       ND     : THE DIMENSION OF THE VECTORS DIF1, DIF2, DIF3, AND ROOT
  59. C
  60. C       N      : THE DEGREE OF THE JACOBI POLYNOMIAL, (i.e. THE NUMBER
  61. C                OF INTERIOR INTERPOLATION POINTS)
  62. C
  63. C       N0     : DETERMINES WHETHER X = 0 IS INCLUDED AS AN
  64. C                INTERPOLATION POINT
  65. C
  66. C                  N0 = 0  ==>  X = 0 IS NOT INCLUDED
  67. C                  N0 = 1  ==>  X = 0 IS INCLUDED
  68. C
  69. C       N1     : DETERMINES WHETHER X = 1 IS INCLUDED AS AN
  70. C                INTERPOLATION POINT
  71. C
  72. C                  N1 = 0  ==>  X = 1 IS NOT INCLUDED
  73. C                  N1 = 1  ==>  X = 1 IS INCLUDED
  74. C
  75. C       ALPHA  : THE VALUE OF ALPHA IN THE DESCRIPTION OF THE JACOBI
  76. C                POLYNOMIAL
  77. C
  78. C       BETA   : THE VALUE OF BETA IN THE DESCRIPTION OF THE JACOBI
  79. C                POLYNOMIAL
  80. C
  81. C       FOR A MORE COMPLETE EXPLANATION OF ALPHA AN BETA, SEE VILLADSEN
  82. C       AND MICHELSEN, PAGES 57 TO 59
  83. C
  84. C     OUTPUT PARAMETERS:
  85. C
  86. C       ROOT   : ONE DIMENSIONAL VECTOR CONTAINING ON EXIT THE
  87. C                N + N0 + N1 ZEROS OF THE NODE POLYNOMIAL USED IN THE
  88. C                INTERPOLATION ROUTINE
  89. C
  90. C       DIF1   : ONE DIMENSIONAL VECTOR CONTAINING THE FIRST DERIVATIVE
  91. C                OF THE NODE POLYNOMIAL AT THE ZEROS
  92. C
  93. C       DIF2   : ONE DIMENSIONAL VECTOR CONTAINING THE SECOND DERIVATIVE
  94. C                OF THE NODE POLYNOMIAL AT THE ZEROS
  95. C
  96. C       DIF3   : ONE DIMENSIONAL VECTOR CONTAINING THE THIRD DERIVATIVE
  97. C                OF THE NODE POLYNOMIAL AT THE ZEROS
  98. C
  99. C     COMMON BLOCKS:      NONE
  100. C
  101. C     REQUIRED ROUTINES:  VILERR, DIF
  102. C
  103. C***********************************************************************
  104. C
  105.       INTEGER           I,J,NT,IER
  106.       DOUBLE PRECISION  AB,AD,AP,Z1,Z,Y,X,XD,XN,XD1,XN1,XP,XP1,ZC
  107.       DOUBLE PRECISION  ZERO,ONE,TWO
  108.       LOGICAL           LSTOP
  109. C
  110.       PARAMETER ( ZERO = 0.0D+00, ONE = 1.0D+00, TWO = 2.0D+00 )
  111. C
  112. C -- ERROR CHECKING
  113. C
  114.       IF ((N0 .NE. 0) .AND. (N0 .NE. 1)) THEN
  115.         IER   = 1
  116.         LSTOP = .TRUE.
  117.         CALL VILERR(IER,LSTOP)
  118.       ELSE
  119.       END IF
  120. C
  121.       IF ((N1 .NE. 0) .AND. (N1 .NE. 1)) THEN
  122.         IER   = 2
  123.         LSTOP = .TRUE.
  124.         CALL VILERR(IER,LSTOP)
  125.       ELSE
  126.       END IF
  127. C
  128.       IF (ND .LT. (N + N0 + N1)) THEN
  129.         IER   = 3
  130.         LSTOP = .TRUE.
  131.         CALL VILERR(IER,LSTOP)
  132.       ELSE
  133.       END IF
  134. C
  135.       IF ((N + N0 + N1) .LT. 1) THEN
  136.         IER   = 7
  137.         LSTOP = .TRUE.
  138.         CALL VILERR(IER,LSTOP)
  139.       ELSE
  140.       END IF
  141. C
  142. C -- FIRST EVALUATION OF COEFFICIENTS IN RECURSION FORMULAS.
  143. C -- RECURSION COEFFICIENTS ARE STORED IN DIF1 AND DIF2.
  144. C
  145.       AB      = ALPHA + BETA
  146.       AD      = BETA - ALPHA
  147.       AP      = BETA*ALPHA
  148.       DIF1(1) = (AD/(AB + TWO) + ONE)/TWO
  149.       DIF2(1) = ZERO
  150. C
  151.       IF(N .GE. 2) THEN
  152.         DO 10 I = 2,N
  153. C
  154.           Z1      = DBLE(I) - ONE
  155.           Z       = AB + 2*Z1
  156.           DIF1(I) = (AB*AD/Z/(Z + TWO) + ONE)/TWO
  157. C
  158.           IF (I .EQ. 2) THEN
  159.             DIF2(I) = (AB + AP + Z1)/Z/Z/(Z + ONE)
  160.           ELSE
  161.             Z       = Z*Z
  162.             Y       = Z1*(AB + Z1)
  163.             Y       = Y*(AP + Y)
  164.             DIF2(I) = Y/Z/(Z - ONE)
  165.           END IF
  166. C
  167.    10   CONTINUE
  168.       ELSE
  169.       END IF
  170. C
  171. C -- ROOT DETERMINATION BY NEWTON METHOD WITH SUPPRESSION OF
  172. C -- PREVIOUSLY DETERMINED ROOTS
  173. C
  174.       X = ZERO
  175. C
  176.       DO 20 I = 1,N
  177. C
  178.    25   CONTINUE
  179.         XD  = ZERO
  180.         XN  = ONE
  181.         XD1 = ZERO
  182.         XN1 = ZERO
  183. C
  184.         DO 30 J = 1,N
  185.           XP  = (DIF1(J) - X)*XN  - DIF2(J)*XD
  186.           XP1 = (DIF1(J) - X)*XN1 - DIF2(J)*XD1 - XN
  187.           XD  = XN
  188.           XD1 = XN1
  189.           XN  = XP
  190.           XN1 = XP1
  191.    30   CONTINUE
  192. C
  193.         ZC  = ONE
  194.         Z   = XN/XN1
  195. C
  196.         IF (I .NE. 1) THEN
  197.           DO 22 J = 2,I
  198.             ZC = ZC - Z/(X - ROOT(J-1))
  199.    22     CONTINUE
  200.         ELSE
  201.         END IF
  202. C
  203.         Z  = Z/ZC
  204.         X  = X - Z
  205. C
  206.         IF (DABS(Z) .GT. 1.D-09) THEN
  207. C
  208. C -- BACKWARD BRANCH
  209. C
  210.           GO TO 25
  211.         ELSE
  212.         END IF
  213. C
  214.         ROOT(I) = X
  215.         X = X + 0.0001D0
  216. C
  217.    20 CONTINUE
  218. C
  219. C -- ADD INTERPOLATION POINTS AT X = 0 AND/OR X = 1
  220. C
  221.       NT = N + N0 + N1
  222. C
  223.       IF (N0 .NE. 0) THEN
  224.         DO 31 I = 1,N
  225.           J = N + 1 - I
  226.           ROOT(J+1) = ROOT(J)
  227.    31   CONTINUE
  228.         ROOT(1) = ZERO
  229.       ELSE
  230.       END IF
  231. C
  232.       IF (N1 .EQ. 1) THEN
  233.         ROOT(NT) = ONE
  234.       ELSE
  235.       END IF
  236. C
  237.       CALL DIF ( NT, ROOT, DIF1, DIF2, DIF3 )
  238. C
  239.       RETURN
  240.       END
  241.