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
|
#include "sys-defines.h"
#include "ode.h"
#include "extern.h"
/*
* Fourth-Order Runge-Kutta
* Copyright Nicholas B. Tufillaro, 1982-1994. All rights reserved.
* GNU enhancements copyright (C) 1996-1997 Free Software Foundation, Inc.
*/
void
#ifdef _HAVE_PROTOS
rk (void)
#else
rk ()
#endif
{
double t;
double halfstep = HALF * tstep;
double onesixth = 1.0 / 6.0;
for (it = 0, t = tstart; !STOPR; t = tstart + (++it) * tstep)
{
symtab->sy_val[0] = symtab->sy_value = t;
field();
for (fsp = dqueue; fsp != NULL; fsp = fsp->sy_link)
{
fsp->sy_val[0] = fsp->sy_value;
fsp->sy_pri[0] = fsp->sy_prime;
}
/* output */
printq();
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]);
}
}
|
