GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
scfode.f
Go to the documentation of this file.
1  SUBROUTINE scfode (METH, ELCO, TESCO)
2 C***BEGIN PROLOGUE SCFODE
3 C***SUBSIDIARY
4 C***PURPOSE Set ODE integrator coefficients.
5 C***TYPE SINGLE PRECISION (SCFODE-S, DCFODE-D)
6 C***AUTHOR Hindmarsh, Alan C., (LLNL)
7 C***DESCRIPTION
8 C
9 C SCFODE is called by the integrator routine to set coefficients
10 C needed there. The coefficients for the current method, as
11 C given by the value of METH, are set for all orders and saved.
12 C The maximum order assumed here is 12 if METH = 1 and 5 if METH = 2.
13 C (A smaller value of the maximum order is also allowed.)
14 C SCFODE is called once at the beginning of the problem,
15 C and is not called again unless and until METH is changed.
16 C
17 C The ELCO array contains the basic method coefficients.
18 C The coefficients el(i), 1 .le. i .le. nq+1, for the method of
19 C order nq are stored in ELCO(i,nq). They are given by a genetrating
20 C polynomial, i.e.,
21 C l(x) = el(1) + el(2)*x + ... + el(nq+1)*x**nq.
22 C For the implicit Adams methods, l(x) is given by
23 C dl/dx = (x+1)*(x+2)*...*(x+nq-1)/factorial(nq-1), l(-1) = 0.
24 C For the BDF methods, l(x) is given by
25 C l(x) = (x+1)*(x+2)* ... *(x+nq)/K,
26 C where K = factorial(nq)*(1 + 1/2 + ... + 1/nq).
27 C
28 C The TESCO array contains test constants used for the
29 C local error test and the selection of step size and/or order.
30 C At order nq, TESCO(k,nq) is used for the selection of step
31 C size at order nq - 1 if k = 1, at order nq if k = 2, and at order
32 C nq + 1 if k = 3.
33 C
34 C***SEE ALSO SLSODE
35 C***ROUTINES CALLED (NONE)
36 C***REVISION HISTORY (YYMMDD)
37 C 791129 DATE WRITTEN
38 C 890501 Modified prologue to SLATEC/LDOC format. (FNF)
39 C 890503 Minor cosmetic changes. (FNF)
40 C 930809 Renamed to allow single/double precision versions. (ACH)
41 C***END PROLOGUE SCFODE
42 C**End
43  INTEGER meth
44  INTEGER i, ib, nq, nqm1, nqp1
45  REAL elco, tesco
46  REAL agamq, fnq, fnqm1, pc, pint, ragq,
47  1 rqfac, rq1fac, tsign, xpin
48  dimension elco(13,12), tesco(3,12)
49  dimension pc(12)
50 C
51 C***FIRST EXECUTABLE STATEMENT SCFODE
52  go to(100, 200), meth
53 C
54  100 elco(1,1) = 1.0e0
55  elco(2,1) = 1.0e0
56  tesco(1,1) = 0.0e0
57  tesco(2,1) = 2.0e0
58  tesco(1,2) = 1.0e0
59  tesco(3,12) = 0.0e0
60  pc(1) = 1.0e0
61  rqfac = 1.0e0
62  DO 140 nq = 2,12
63 C-----------------------------------------------------------------------
64 C The PC array will contain the coefficients of the polynomial
65 C p(x) = (x+1)*(x+2)*...*(x+nq-1).
66 C Initially, p(x) = 1.
67 C-----------------------------------------------------------------------
68  rq1fac = rqfac
69  rqfac = rqfac/nq
70  nqm1 = nq - 1
71  fnqm1 = nqm1
72  nqp1 = nq + 1
73 C Form coefficients of p(x)*(x+nq-1). ----------------------------------
74  pc(nq) = 0.0e0
75  DO 110 ib = 1,nqm1
76  i = nqp1 - ib
77  110 pc(i) = pc(i-1) + fnqm1*pc(i)
78  pc(1) = fnqm1*pc(1)
79 C Compute integral, -1 to 0, of p(x) and x*p(x). -----------------------
80  pint = pc(1)
81  xpin = pc(1)/2.0e0
82  tsign = 1.0e0
83  DO 120 i = 2,nq
84  tsign = -tsign
85  pint = pint + tsign*pc(i)/i
86  120 xpin = xpin + tsign*pc(i)/(i+1)
87 C Store coefficients in ELCO and TESCO. --------------------------------
88  elco(1,nq) = pint*rq1fac
89  elco(2,nq) = 1.0e0
90  DO 130 i = 2,nq
91  130 elco(i+1,nq) = rq1fac*pc(i)/i
92  agamq = rqfac*xpin
93  ragq = 1.0e0/agamq
94  tesco(2,nq) = ragq
95  IF (nq .LT. 12) tesco(1,nqp1) = ragq*rqfac/nqp1
96  tesco(3,nqm1) = ragq
97  140 CONTINUE
98  RETURN
99 C
100  200 pc(1) = 1.0e0
101  rq1fac = 1.0e0
102  DO 230 nq = 1,5
103 C-----------------------------------------------------------------------
104 C The PC array will contain the coefficients of the polynomial
105 C p(x) = (x+1)*(x+2)*...*(x+nq).
106 C Initially, p(x) = 1.
107 C-----------------------------------------------------------------------
108  fnq = nq
109  nqp1 = nq + 1
110 C Form coefficients of p(x)*(x+nq). ------------------------------------
111  pc(nqp1) = 0.0e0
112  DO 210 ib = 1,nq
113  i = nq + 2 - ib
114  210 pc(i) = pc(i-1) + fnq*pc(i)
115  pc(1) = fnq*pc(1)
116 C Store coefficients in ELCO and TESCO. --------------------------------
117  DO 220 i = 1,nqp1
118  220 elco(i,nq) = pc(i)/pc(2)
119  elco(2,nq) = 1.0e0
120  tesco(1,nq) = rq1fac
121  tesco(2,nq) = nqp1/elco(1,nq)
122  tesco(3,nq) = (nq+2)/elco(1,nq)
123  rq1fac = rq1fac/fnq
124  230 CONTINUE
125  RETURN
126 C----------------------- END OF SUBROUTINE SCFODE ----------------------
127  END