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
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