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
113
114
115
116
117
118
119
120
121
122
|
#include "sys-defines.h"
#include "ode.h"
#include "extern.h"
#define PASTVAL (3) /* previous values, val[0] is current value */
/*
* Adams-Moulton with constant step size
* Copyright Nicholas B. Tufillaro, 1982-1994. All rights reserved.
* GNU enhancements copyright (C) 1996-1997 Free Software Foundation, Inc.
*/
void
#ifdef _HAVE_PROTOS
am (void)
#else
am ()
#endif
{
double t;
double halfstep = HALF * tstep;
double sconst = tstep / 24.0; /* step constant */
double onesixth = 1.0 / 6.0;
/* Runge-Kutta startup */
for (it = 0, t = tstart; it <= PASTVAL && !STOPR; t = tstart + (++it) * tstep)
{
symtab->sy_value = symtab->sy_val[0] = t;
field();
for (fsp = dqueue; fsp != NULL; fsp = fsp->sy_link)
{
int j;
for (j = it; j > 0; j--)
{
fsp->sy_val[j] = fsp->sy_val[j-1];
fsp->sy_pri[j] = fsp->sy_pri[j-1];
}
fsp->sy_pri[0] = fsp->sy_prime;
fsp->sy_val[0] = fsp->sy_value;
}
/* output */
printq();
if (it == PASTVAL)
break; /* startup complete */
for (fsp = dqueue; fsp != NULL; fsp = fsp->sy_link)
{
fsp->sy_k[0] = tstep * fsp->sy_prime;
fsp->sy_value = fsp->sy_val[0] + HALF * fsp->sy_k[0];
}
symtab->sy_value = t + halfstep;
field();
for (fsp = dqueue; fsp != NULL; fsp = fsp->sy_link)
{
fsp->sy_k[1] = tstep * fsp->sy_prime;
fsp->sy_value = fsp->sy_val[0] + HALF * fsp->sy_k[1];
}
symtab->sy_value = t + halfstep;
field();
for (fsp = dqueue; fsp != NULL; fsp = fsp->sy_link)
{
fsp->sy_k[2] = tstep * fsp->sy_prime;
fsp->sy_value = fsp->sy_val[0] + fsp->sy_k[2];
}
symtab->sy_value = t + tstep;
field();
for (fsp = dqueue; fsp != NULL; fsp = fsp->sy_link)
fsp->sy_k[3] = tstep * fsp->sy_prime;
for (fsp = dqueue; fsp != NULL; fsp = fsp->sy_link)
{
fsp->sy_value = fsp->sy_val[0]
+ onesixth * (fsp->sy_k[0]
+ TWO * fsp->sy_k[1]
+ TWO * fsp->sy_k[2]
+ fsp->sy_k[3]);
}
}
/* predictor - corrector */
while (!STOPA)
{
/* Adams-Bashforth predictor */
for (fsp = dqueue; fsp != NULL ; fsp = fsp->sy_link)
{
fsp->sy_value = fsp->sy_val[0]
+ (sconst) * (55 * fsp->sy_pri[0]
-59 * fsp->sy_pri[1]
+37 * fsp->sy_pri[2]
-9 * fsp->sy_pri[3]);
}
symtab->sy_val[0] = symtab->sy_value = t = tstart + (++it) * tstep;
field();
/* Adams-Moulton corrector */
for (fsp = dqueue; fsp != NULL; fsp = fsp->sy_link)
{
fsp->sy_value = fsp->sy_val[0]
+ (sconst) * (9 * fsp->sy_prime
+19 * fsp->sy_pri[0]
-5 * fsp->sy_pri[1]
+ fsp->sy_pri[2]);
}
field();
/* cycle indices */
for (fsp = dqueue; fsp != NULL; fsp = fsp->sy_link)
{
int j;
for (j = PASTVAL; j > 0; j--)
{
fsp->sy_val[j] = fsp->sy_val[j-1];
fsp->sy_pri[j] = fsp->sy_pri[j-1];
}
fsp->sy_val[0] = fsp->sy_value;
fsp->sy_pri[0] = fsp->sy_prime;
}
/* output */
printq();
}
}
|
