cfode.f

Go to the documentation of this file.

      SUBROUTINE cfode (METH, ELCO, TESCO)
CLLL. OPTIMIZE
      INTEGER meth
      INTEGER i, ib, nq, nqm1, nqp1
      DOUBLE PRECISION elco, tesco
      DOUBLE PRECISION agamq, fnq, fnqm1, pc, pint, ragq,
     1   rqfac, rq1fac, tsign, xpin
      dimension elco(13,12), tesco(3,12)
C-----------------------------------------------------------------------
C CFODE IS CALLED BY THE INTEGRATOR ROUTINE TO SET COEFFICIENTS
C NEEDED THERE.  THE COEFFICIENTS FOR THE CURRENT METHOD, AS
C GIVEN BY THE VALUE OF METH, ARE SET FOR ALL ORDERS AND SAVED.
C THE MAXIMUM ORDER ASSUMED HERE IS 12 IF METH = 1 AND 5 IF METH = 2. 
C (A SMALLER VALUE OF THE MAXIMUM ORDER IS ALSO ALLOWED.)
C CFODE IS CALLED ONCE AT THE BEGINNING OF THE PROBLEM,
C AND IS NOT CALLED AGAIN UNLESS AND UNTIL METH IS CHANGED. 
C
C THE ELCO ARRAY CONTAINS THE BASIC METHOD COEFFICIENTS.
C THE COEFFICIENTS EL(I), 1 .LE. I .LE. NQ+1, FOR THE METHOD OF
C ORDER NQ ARE STORED IN ELCO(I,NQ).  THEY ARE GIVEN BY A GENETRATING 
C POLYNOMIAL, I.E., 
C     L(X) = EL(1) + EL(2)*X + ... + EL(NQ+1)*X**NQ.
C FOR THE IMPLICIT ADAMS METHODS, L(X) IS GIVEN BY
C     DL/DX = (X+1)*(X+2)*...*(X+NQ-1)/FACTORIAL(NQ-1),    L(-1) = 0. 
C FOR THE BDF METHODS, L(X) IS GIVEN BY 
C     L(X) = (X+1)*(X+2)* ... *(X+NQ)/K,
C WHERE         K = FACTORIAL(NQ)*(1 + 1/2 + ... + 1/NQ).
C
C THE TESCO ARRAY CONTAINS TEST CONSTANTS USED FOR THE
C LOCAL ERROR TEST AND THE SELECTION OF STEP SIZE AND/OR ORDER.
C AT ORDER NQ, TESCO(K,NQ) IS USED FOR THE SELECTION OF STEP
C SIZE AT ORDER NQ - 1 IF K = 1, AT ORDER NQ IF K = 2, AND AT ORDER
C NQ + 1 IF K = 3.
C-----------------------------------------------------------------------
      dimension pc(12)
C
      go to(100, 200), meth
C
 100  elco(1,1) = 1.0d0
      elco(2,1) = 1.0d0
      tesco(1,1) = 0.0d0
      tesco(2,1) = 2.0d0
      tesco(1,2) = 1.0d0
      tesco(3,12) = 0.0d0
      pc(1) = 1.0d0 
      rqfac = 1.0d0 
      DO 140 nq = 2,12
C-----------------------------------------------------------------------
C THE PC ARRAY WILL CONTAIN THE COEFFICIENTS OF THE POLYNOMIAL
C     P(X) = (X+1)*(X+2)*...*(X+NQ-1).
C INITIALLY, P(X) = 1.
C-----------------------------------------------------------------------
        rq1fac = rqfac
        rqfac = rqfac/dble(nq)
        nqm1 = nq - 1
        fnqm1 = dble(nqm1)  
        nqp1 = nq + 1
C FORM COEFFICIENTS OF P(X)*(X+NQ-1). ----------------------------------
        pc(nq) = 0.0d0
        DO 110 ib = 1,nqm1
          i = nqp1 - ib
 110      pc(i) = pc(i-1) + fnqm1*pc(i) 
        pc(1) = fnqm1*pc(1)
C COMPUTE INTEGRAL, -1 TO 0, OF P(X) AND X*P(X). -----------------------
        pint = pc(1)
        xpin = pc(1)/2.0d0
        tsign = 1.0d0
        DO 120 i = 2,nq
          tsign = -tsign
          pint = pint + tsign*pc(i)/dble(i)     
 120      xpin = xpin + tsign*pc(i)/dble(i+1)   
C STORE COEFFICIENTS IN ELCO AND TESCO. --------------------------------
        elco(1,nq) = pint*rq1fac
        elco(2,nq) = 1.0d0
        DO 130 i = 2,nq
 130      elco(i+1,nq) = rq1fac*pc(i)/dble(i)   
        agamq = rqfac*xpin
        ragq = 1.0d0/agamq
        tesco(2,nq) = ragq
        IF (nq .LT. 12) tesco(1,nqp1) = ragq*rqfac/dble(nqp1)       
        tesco(3,nqm1) = ragq
 140    CONTINUE
      RETURN
C
 200  pc(1) = 1.0d0 
      rq1fac = 1.0d0
      DO 230 nq = 1,5
C-----------------------------------------------------------------------
C THE PC ARRAY WILL CONTAIN THE COEFFICIENTS OF THE POLYNOMIAL
C     P(X) = (X+1)*(X+2)*...*(X+NQ).
C INITIALLY, P(X) = 1.
C-----------------------------------------------------------------------
        fnq = dble(nq)      
        nqp1 = nq + 1
C FORM COEFFICIENTS OF P(X)*(X+NQ). ------------------------------------
        pc(nqp1) = 0.0d0
        DO 210 ib = 1,nq
          i = nq + 2 - ib
 210      pc(i) = pc(i-1) + fnq*pc(i)
        pc(1) = fnq*pc(1)
C STORE COEFFICIENTS IN ELCO AND TESCO. --------------------------------
        DO 220 i = 1,nqp1
 220      elco(i,nq) = pc(i)/pc(2)
        elco(2,nq) = 1.0d0
        tesco(1,nq) = rq1fac
        tesco(2,nq) = dble(nqp1)/elco(1,nq)     
        tesco(3,nq) = dble(nq+2)/elco(1,nq)     
        rq1fac = rq1fac/fnq
 230    CONTINUE
      RETURN
C----------------------- END OF SUBROUTINE CFODE -----------------------
      END