GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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
scfode
subroutine scfode(METH, ELCO, TESCO)
Definition: scfode.f:2