GNU Octave  4.2.1 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
cfode.f
Go to the documentation of this file.
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
may be zero for pure relative error test tem the relative tolerance must be greater than or equal to