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 / odepack / cfode.f next >
Text File  |  1996-09-28  |  4KB  |  113 lines
  1.       SUBROUTINE CFODE (METH, ELCO, TESCO)
  2. CLLL. OPTIMIZE
  3.       INTEGER METH
  4.       INTEGER I, IB, NQ, NQM1, NQP1
  5.       DOUBLE PRECISION ELCO, TESCO
  6.       DOUBLE PRECISION AGAMQ, FNQ, FNQM1, PC, PINT, RAGQ,
  7.      1   RQFAC, RQ1FAC, TSIGN, XPIN
  8.       DIMENSION ELCO(13,12), TESCO(3,12)
  9. C-----------------------------------------------------------------------
  10. C CFODE IS CALLED BY THE INTEGRATOR ROUTINE TO SET COEFFICIENTS
  11. C NEEDED THERE.  THE COEFFICIENTS FOR THE CURRENT METHOD, AS
  12. C GIVEN BY THE VALUE OF METH, ARE SET FOR ALL ORDERS AND SAVED.
  13. C THE MAXIMUM ORDER ASSUMED HERE IS 12 IF METH = 1 AND 5 IF METH = 2. 
  14. C (A SMALLER VALUE OF THE MAXIMUM ORDER IS ALSO ALLOWED.)
  15. C CFODE IS CALLED ONCE AT THE BEGINNING OF THE PROBLEM,
  16. C AND IS NOT CALLED AGAIN UNLESS AND UNTIL METH IS CHANGED. 
  17. C
  18. C THE ELCO ARRAY CONTAINS THE BASIC METHOD COEFFICIENTS.
  19. C THE COEFFICIENTS EL(I), 1 .LE. I .LE. NQ+1, FOR THE METHOD OF
  20. C ORDER NQ ARE STORED IN ELCO(I,NQ).  THEY ARE GIVEN BY A GENETRATING 
  21. C POLYNOMIAL, I.E., 
  22. C     L(X) = EL(1) + EL(2)*X + ... + EL(NQ+1)*X**NQ.
  23. C FOR THE IMPLICIT ADAMS METHODS, L(X) IS GIVEN BY
  24. C     DL/DX = (X+1)*(X+2)*...*(X+NQ-1)/FACTORIAL(NQ-1),    L(-1) = 0. 
  25. C FOR THE BDF METHODS, L(X) IS GIVEN BY 
  26. C     L(X) = (X+1)*(X+2)* ... *(X+NQ)/K,
  27. C WHERE         K = FACTORIAL(NQ)*(1 + 1/2 + ... + 1/NQ).
  28. C
  29. C THE TESCO ARRAY CONTAINS TEST CONSTANTS USED FOR THE
  30. C LOCAL ERROR TEST AND THE SELECTION OF STEP SIZE AND/OR ORDER.
  31. C AT ORDER NQ, TESCO(K,NQ) IS USED FOR THE SELECTION OF STEP
  32. C SIZE AT ORDER NQ - 1 IF K = 1, AT ORDER NQ IF K = 2, AND AT ORDER
  33. C NQ + 1 IF K = 3.
  34. C-----------------------------------------------------------------------
  35.       DIMENSION PC(12)
  36. C
  37.       GO TO (100, 200), METH
  38. C
  39.  100  ELCO(1,1) = 1.0D0
  40.       ELCO(2,1) = 1.0D0
  41.       TESCO(1,1) = 0.0D0
  42.       TESCO(2,1) = 2.0D0
  43.       TESCO(1,2) = 1.0D0
  44.       TESCO(3,12) = 0.0D0
  45.       PC(1) = 1.0D0 
  46.       RQFAC = 1.0D0 
  47.       DO 140 NQ = 2,12
  48. C-----------------------------------------------------------------------
  49. C THE PC ARRAY WILL CONTAIN THE COEFFICIENTS OF THE POLYNOMIAL
  50. C     P(X) = (X+1)*(X+2)*...*(X+NQ-1).
  51. C INITIALLY, P(X) = 1.
  52. C-----------------------------------------------------------------------
  53.         RQ1FAC = RQFAC
  54.         RQFAC = RQFAC/DBLE(NQ)
  55.         NQM1 = NQ - 1
  56.         FNQM1 = DBLE(NQM1)  
  57.         NQP1 = NQ + 1
  58. C FORM COEFFICIENTS OF P(X)*(X+NQ-1). ----------------------------------
  59.         PC(NQ) = 0.0D0
  60.         DO 110 IB = 1,NQM1
  61.           I = NQP1 - IB
  62.  110      PC(I) = PC(I-1) + FNQM1*PC(I) 
  63.         PC(1) = FNQM1*PC(1)
  64. C COMPUTE INTEGRAL, -1 TO 0, OF P(X) AND X*P(X). -----------------------
  65.         PINT = PC(1)
  66.         XPIN = PC(1)/2.0D0
  67.         TSIGN = 1.0D0
  68.         DO 120 I = 2,NQ
  69.           TSIGN = -TSIGN
  70.           PINT = PINT + TSIGN*PC(I)/DBLE(I)     
  71.  120      XPIN = XPIN + TSIGN*PC(I)/DBLE(I+1)   
  72. C STORE COEFFICIENTS IN ELCO AND TESCO. --------------------------------
  73.         ELCO(1,NQ) = PINT*RQ1FAC
  74.         ELCO(2,NQ) = 1.0D0
  75.         DO 130 I = 2,NQ
  76.  130      ELCO(I+1,NQ) = RQ1FAC*PC(I)/DBLE(I)   
  77.         AGAMQ = RQFAC*XPIN
  78.         RAGQ = 1.0D0/AGAMQ
  79.         TESCO(2,NQ) = RAGQ
  80.         IF (NQ .LT. 12) TESCO(1,NQP1) = RAGQ*RQFAC/DBLE(NQP1)       
  81.         TESCO(3,NQM1) = RAGQ
  82.  140    CONTINUE
  83.       RETURN
  84. C
  85.  200  PC(1) = 1.0D0 
  86.       RQ1FAC = 1.0D0
  87.       DO 230 NQ = 1,5
  88. C-----------------------------------------------------------------------
  89. C THE PC ARRAY WILL CONTAIN THE COEFFICIENTS OF THE POLYNOMIAL
  90. C     P(X) = (X+1)*(X+2)*...*(X+NQ).
  91. C INITIALLY, P(X) = 1.
  92. C-----------------------------------------------------------------------
  93.         FNQ = DBLE(NQ)      
  94.         NQP1 = NQ + 1
  95. C FORM COEFFICIENTS OF P(X)*(X+NQ). ------------------------------------
  96.         PC(NQP1) = 0.0D0
  97.         DO 210 IB = 1,NQ
  98.           I = NQ + 2 - IB
  99.  210      PC(I) = PC(I-1) + FNQ*PC(I)
  100.         PC(1) = FNQ*PC(1)
  101. C STORE COEFFICIENTS IN ELCO AND TESCO. --------------------------------
  102.         DO 220 I = 1,NQP1
  103.  220      ELCO(I,NQ) = PC(I)/PC(2)
  104.         ELCO(2,NQ) = 1.0D0
  105.         TESCO(1,NQ) = RQ1FAC
  106.         TESCO(2,NQ) = DBLE(NQP1)/ELCO(1,NQ)     
  107.         TESCO(3,NQ) = DBLE(NQ+2)/ELCO(1,NQ)     
  108.         RQ1FAC = RQ1FAC/FNQ
  109.  230    CONTINUE
  110.       RETURN
  111. C----------------------- END OF SUBROUTINE CFODE -----------------------
  112.       END 
  113.