1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
|
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/dfloat(nq)
nqm1 = nq - 1
fnqm1 = dfloat(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)/dfloat(i)
120 xpin = xpin + tsign*pc(i)/dfloat(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)/dfloat(i)
agamq = rqfac*xpin
ragq = 1.0d0/agamq
tesco(2,nq) = ragq
if (nq .lt. 12) tesco(1,nqp1) = ragq*rqfac/dfloat(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 = dfloat(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) = dfloat(nqp1)/elco(1,nq)
tesco(3,nq) = dfloat(nq+2)/elco(1,nq)
rq1fac = rq1fac/fnq
230 continue
return
c----------------------- end of subroutine cfode -----------------------
end
|
