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 / radau.f < prev    next >
Text File  |  1996-09-28  |  5KB  |  210 lines
  1.       SUBROUTINE RADAU
  2.      +  (
  3.      +  ND, N, N0, N1, ID, ALPHA, BETA, ROOT, DIF1, VECT
  4.      +  )
  5. C
  6.       INTEGER           ND, N, N0, N1, ID
  7.       DOUBLE PRECISION  ALPHA, BETA, ROOT(ND), DIF1(ND), VECT(ND)
  8. C
  9. C***********************************************************************
  10. C
  11. C     RADAU OR LOBATTO QUADRATURE
  12. C
  13. C     VILLADSEN AND MICHELSEN, PAGES 133-135, 419
  14. C
  15. C     INPUT PARAMETERS:
  16. C
  17. C       ND     : THE DIMENSION OF THE VECTORS DIF1, DIF2, DIF3, AND ROOT
  18. C
  19. C       N      : THE DEGREE OF THE JACOBI POLYNOMIAL, (i.e. THE NUMBER
  20. C                OF INTERIOR INTERPOLATION POINTS)
  21. C
  22. C       N0     : DETERMINES WHETHER X = 0 IS INCLUDED AS AN
  23. C                INTERPOLATION POINT
  24. C
  25. C                  N0 = 0  ==>  X = 0 IS NOT INCLUDED
  26. C                  N0 = 1  ==>  X = 0 IS INCLUDED
  27. C
  28. C       N1     : DETERMINES WHETHER X = 1 IS INCLUDED AS AN
  29. C                INTERPOLATION POINT
  30. C
  31. C                  N1 = 0  ==>  X = 1 IS NOT INCLUDED
  32. C                  N1 = 1  ==>  X = 1 IS INCLUDED
  33. C
  34. C       ID     : INDICATOR
  35. C
  36. C                  ID = 1  ==>  RADAU QUADRATURE WEIGHTS INCLUDING X = 1
  37. C                  ID = 2  ==>  RADAU QUADRATURE WEIGHTS INCLUDING X = 0
  38. C                  ID = 3  ==>  LOBATTO QUADRATURE WEIGHTS INCLUDING
  39. C                               BOTH X = 0 AND X = 1
  40. C
  41. C       ALPHA  : THE VALUE OF ALPHA IN THE DESCRIPTION OF THE JACOBI
  42. C                POLYNOMIAL
  43. C
  44. C       BETA   : THE VALUE OF BETA IN THE DESCRIPTION OF THE JACOBI
  45. C                POLYNOMIAL
  46. C
  47. C                FOR A MORE COMPLETE EXPLANATION OF ALPHA AN BETA, SEE
  48. C                VILLADSEN AND MICHELSEN, PAGES 57 TO 59
  49. C
  50. C       ROOT   : ONE DIMENSIONAL VECTOR CONTAINING ON EXIT THE
  51. C                N + N0 + N1 ZEROS OF THE NODE POLYNOMIAL USED IN THE
  52. C                INTERPOLATION ROUTINE
  53. C
  54. C       DIF1   : ONE DIMENSIONAL VECTOR CONTAINING THE FIRST DERIVATIVE
  55. C                OF THE NODE POLYNOMIAL AT THE ZEROS
  56. C
  57. C       THE NODE POLYNOMIAL IS GIVEN BY
  58. C
  59. C                     N0    (ALPHA',BETA')          N1
  60. C         P  (X)  =  X   * P (X)           * (X - 1)
  61. C          NT               N
  62. C
  63. C       THE ARGUMENTS ALPHA' AND BETA' TO BE USED IN JCOBI FOR
  64. C       CALCULATION OF ROOT AND DIF1 DEPEND ON WHETHER X = 0 , X = 1 OR
  65. C       BOTH ARE USED AS EXTRA QUADRATURE POINTS.  THUS:
  66. C
  67. C         ID = 1:  ALPHA' = ALPHA + 1, BETA' = BETA
  68. C         ID = 2:  ALPHA' = ALPHA    , BETA' = BETA + 1
  69. C         ID = 3:  ALPHA' = ALPHA + 1, BETA' = BETA + 1
  70. C
  71. C       NOTE:
  72. C
  73. C         ID = 1  REQUIRES THAT N0 = 0 OR 1, N1 = 1
  74. C         ID = 2  REQUIRES THAT N0 = 1     , N1 = 0 OR 1
  75. C         ID = 3  REQUIRES THAT N0 = 1     , N1 = 1
  76. C
  77. C     OUTPUT PARAMETERS:
  78. C
  79. C       VECT   : VECTOR OF THE NT COMPUTED QUADRATURE WEIGHTS,
  80. C                NORMALIZED SUCH THAT
  81. C
  82. C                   NT
  83. C                  SUM  VECT(I) = 1
  84. C                  I=1
  85. C
  86. C                FOR A MORE COMPLETE EXPLANATION SEE VILLADSEN AND
  87. C                MICHELSEN, PAGES 133 TO 135
  88. C
  89. C     COMMON BLOCKS:      NONE
  90. C
  91. C     REQUIRED ROUTINES:  VILERR
  92. C
  93. C***********************************************************************
  94. C
  95.       INTEGER           I,NT,IER
  96.       DOUBLE PRECISION  AX,S,X
  97.       DOUBLE PRECISION  ZERO,ONE
  98.       LOGICAL           LSTOP
  99. C
  100.       PARAMETER ( ZERO = 0.0D+00, ONE = 1.0D+00 )
  101. C
  102. C -- ERROR CHECKING
  103. C
  104.       IF ((N0 .NE. 0) .AND. (N0 .NE. 1)) THEN
  105.         IER   = 1
  106.         LSTOP = .TRUE.
  107.         CALL VILERR(IER,LSTOP)
  108.       ELSE
  109.       END IF
  110. C
  111.       IF ((N1 .NE. 0) .AND. (N1 .NE. 1)) THEN
  112.         IER   = 2
  113.         LSTOP = .TRUE.
  114.         CALL VILERR(IER,LSTOP)
  115.       ELSE
  116.       END IF
  117. C
  118.       IF (ND .LT. (N + N0 + N1)) THEN
  119.         IER   = 3
  120.         LSTOP = .TRUE.
  121.         CALL VILERR(IER,LSTOP)
  122.       ELSE
  123.       END IF
  124. C
  125.       IF ((N + N0 + N1) .LT. 1) THEN
  126.         IER   = 7
  127.         LSTOP = .TRUE.
  128.         CALL VILERR(IER,LSTOP)
  129.       ELSE
  130.       END IF
  131. C
  132.       IF ((ID .NE. 1) .AND. (ID.NE. 2) .AND. (ID .NE. 3)) THEN
  133.         IER   = 8
  134.         LSTOP = .TRUE.
  135.         CALL VILERR(IER,LSTOP)
  136.       ELSE
  137.       END IF
  138. C
  139.       IF ((ID .EQ. 1) .AND. (N1 .NE. 1)) THEN
  140.         IER   = 9
  141.         LSTOP = .TRUE.
  142.         CALL VILERR(IER,LSTOP)
  143.       ELSE
  144.       END IF
  145. C
  146.       IF ((ID .EQ. 2) .AND. (N0 .NE. 1)) THEN
  147.         IER   = 10
  148.         LSTOP = .TRUE.
  149.         CALL VILERR(IER,LSTOP)
  150.       ELSE
  151.       END IF
  152. C
  153.       IF ((ID .EQ. 3) .AND. ((N0 .NE. 1) .OR. (N1 .NE. 1))) THEN
  154.         IER   = 11
  155.         LSTOP = .TRUE.
  156.         CALL VILERR(IER,LSTOP)
  157.       ELSE
  158.       END IF
  159. C
  160. C -- EVALUATE RADAU OR LOBATTO QUADRATURE WEIGHTS
  161. C
  162.       S  = ZERO
  163.       NT = N + N0 + N1
  164. C
  165.       DO 40 I = 1,NT
  166. C
  167.         X = ROOT(I)
  168. C
  169.         IF      (ID .EQ. 1) THEN
  170.           AX = X
  171.           IF (N0 .EQ. 0) THEN
  172.             AX = ONE/AX
  173.           ELSE
  174.           END IF
  175.         ELSE IF (ID .EQ. 2) THEN
  176.           AX = ONE - X
  177.           IF (N1 .EQ. 0) THEN
  178.             AX = ONE/AX
  179.           ELSE
  180.           END IF
  181.         ELSE IF (ID .EQ. 3) THEN
  182.           AX = ONE
  183.         ELSE
  184.         END IF
  185. C
  186.         VECT(I) = AX/DIF1(I)**2
  187. C
  188.    40 CONTINUE
  189. C
  190.       IF (ID .NE. 2) THEN
  191.         VECT(NT) = VECT(NT)/(ONE + ALPHA)
  192.       ELSE
  193.       END IF
  194. C
  195.       IF (ID .GT. 1) THEN
  196.         VECT(1)  = VECT( 1)/(ONE + BETA)
  197.       ELSE
  198.       END IF
  199. C
  200.       DO 50 I = 1,NT
  201.         S = S + VECT(I)
  202.    50 CONTINUE
  203. C
  204.       DO 60 I = 1,NT
  205.         VECT(I) = VECT(I)/S
  206.    60 CONTINUE
  207. C
  208.       RETURN
  209.       END
  210.