line |
stmt |
bran |
cond |
sub |
pod |
time |
code |
1
|
|
|
|
|
|
|
package Text::NumericData::App::txdodeint; |
2
|
|
|
|
|
|
|
|
3
|
1
|
|
|
1
|
|
346
|
use Text::NumericData::App; |
|
1
|
|
|
|
|
2
|
|
|
1
|
|
|
|
|
27
|
|
4
|
1
|
|
|
1
|
|
5
|
use Text::NumericData::Calc qw(formula_function); |
|
1
|
|
|
|
|
2
|
|
|
1
|
|
|
|
|
37
|
|
5
|
|
|
|
|
|
|
|
6
|
1
|
|
|
1
|
|
4
|
use strict; |
|
1
|
|
|
|
|
1
|
|
|
1
|
|
|
|
|
1289
|
|
7
|
|
|
|
|
|
|
|
8
|
|
|
|
|
|
|
# This is just a placeholder because of a past build system bug. |
9
|
|
|
|
|
|
|
# The one and only version for Text::NumericData is kept in |
10
|
|
|
|
|
|
|
# the Text::NumericData module itself. |
11
|
|
|
|
|
|
|
our $VERSION = '1'; |
12
|
|
|
|
|
|
|
$VERSION = eval $VERSION; |
13
|
|
|
|
|
|
|
|
14
|
|
|
|
|
|
|
my $infostring = "integrate a given ordinary differential equation system along a coordinate |
15
|
|
|
|
|
|
|
|
16
|
|
|
|
|
|
|
Usage: |
17
|
|
|
|
|
|
|
pipe | txdodeint [parameters] [ [ [initval ...]]] | pipe |
18
|
|
|
|
|
|
|
|
19
|
|
|
|
|
|
|
The first fixed parameter is the column to use as time coordinate for which |
20
|
|
|
|
|
|
|
the expression(s) of the ODE are defined. The derivative can be a |
21
|
|
|
|
|
|
|
system of ODEs. The time derivative values shall be stored in the array members |
22
|
|
|
|
|
|
|
D0 .. Dn, with the current corresponding variable values available in V0 .. Vn |
23
|
|
|
|
|
|
|
and the corresponding values of the (interpolated) auxilliary timeseries |
24
|
|
|
|
|
|
|
from the piped data in [1] .. [m]. There are also the auxillary array values |
25
|
|
|
|
|
|
|
A0 .. Ak (to be used at your leasure to store/load values) and the constants |
26
|
|
|
|
|
|
|
C0 .. Cl (initialised from program parameters). |
27
|
|
|
|
|
|
|
|
28
|
|
|
|
|
|
|
For convenience, the basic setup of time column, ODE, and intial values can |
29
|
|
|
|
|
|
|
be given directly on the command line without mentioning parameter names. |
30
|
|
|
|
|
|
|
The integrated values are written out at the points in time given by the input |
31
|
|
|
|
|
|
|
data. For each variable (size of the array V), a column is appended. |
32
|
|
|
|
|
|
|
|
33
|
|
|
|
|
|
|
Example for a simple system (constant acceleration): |
34
|
|
|
|
|
|
|
|
35
|
|
|
|
|
|
|
D0=V1; D1=A0; D2=C0*[1] |
36
|
|
|
|
|
|
|
|
37
|
|
|
|
|
|
|
Assuming the time in column 1 and the constant acceleration in C0, |
38
|
|
|
|
|
|
|
this computes the evolution of the covered distance V0 via the numerically |
39
|
|
|
|
|
|
|
accelerated speed V1 and also directly from the analytically accelerated |
40
|
|
|
|
|
|
|
speed as variable V2, to give a hint about the accuracy of the numerical |
41
|
|
|
|
|
|
|
integration. |
42
|
|
|
|
|
|
|
|
43
|
|
|
|
|
|
|
The employed integration method is a standard Runge-Kutta scheme with up to |
44
|
|
|
|
|
|
|
4 stages, which should be fine for any application where you consider a |
45
|
|
|
|
|
|
|
humble Perl script for your numerical integration. A comparison of the fully |
46
|
|
|
|
|
|
|
numerical with the fully analytical solution can be constructed via the |
47
|
|
|
|
|
|
|
pipeline |
48
|
|
|
|
|
|
|
|
49
|
|
|
|
|
|
|
txdconstruct -n=11 '[1]=C0-1' \\ |
50
|
|
|
|
|
|
|
| txdodeint --rksteps=1 'D0=V1; D1=3' 0 0 \\ |
51
|
|
|
|
|
|
|
| txdcalc '[4]=3/2*[1]**2; [5]=[4] ? ([2]-[4])/[4] : 0' \\ |
52
|
|
|
|
|
|
|
| txdfilter -N=%g 'integration test' 't/s' 's/m' 'v/m' 's_ref/m' 'error' |
53
|
|
|
|
|
|
|
|
54
|
|
|
|
|
|
|
With rksteps>1, you will not see any difference in this example. In general, |
55
|
|
|
|
|
|
|
the choice of integration stages, time step and interpolation might have |
56
|
|
|
|
|
|
|
significant influence on your results. A simple test of the quality of the |
57
|
|
|
|
|
|
|
chosen integration employs trivial polynomials: |
58
|
|
|
|
|
|
|
|
59
|
|
|
|
|
|
|
txdconstruct -n=11 '[1]=C0-1' \\ |
60
|
|
|
|
|
|
|
| txdodeint --rksteps=2 --timediv=1 \\ |
61
|
|
|
|
|
|
|
'D0=1; D1=2*[1]; D2=3*[1]**2; D3=4*[1]**3' \\ |
62
|
|
|
|
|
|
|
0 0 0 0 \\ |
63
|
|
|
|
|
|
|
| txdcalc '[2]-=[1]; [3]-=[1]**2; [4]-=[1]**3; [5]-=[1]**4' \\ |
64
|
|
|
|
|
|
|
| txdfilter -N=g 'integration order test' x err0 err1 err2 err3 |
65
|
|
|
|
|
|
|
|
66
|
|
|
|
|
|
|
These polynomials can actually be solved exactly to machine precision |
67
|
|
|
|
|
|
|
(depending on rksteps value) and smaller time steps would introduce |
68
|
|
|
|
|
|
|
rounding errors here from the summation. Finally, another classic example, |
69
|
|
|
|
|
|
|
the Lorenz attractor: |
70
|
|
|
|
|
|
|
|
71
|
|
|
|
|
|
|
txdconstruct -n=5001 '[1]=(C0-1)/100' | txdodeint --timediv=1 \\ |
72
|
|
|
|
|
|
|
'D0=10*(V1-V0); D1=28*V0-V1-V0*V2; D2=-8/3*V2+V0*V1' \\ |
73
|
|
|
|
|
|
|
20 -20 1 |
74
|
|
|
|
|
|
|
|
75
|
|
|
|
|
|
|
More practical applications actually have some more data columns besides |
76
|
|
|
|
|
|
|
the time in the input data (measurements) and involve derivative expressions |
77
|
|
|
|
|
|
|
that make use of this data-driven time-dependence."; |
78
|
|
|
|
|
|
|
|
79
|
|
|
|
|
|
|
our @ISA = ('Text::NumericData::App'); |
80
|
|
|
|
|
|
|
|
81
|
|
|
|
|
|
|
sub new |
82
|
|
|
|
|
|
|
{ |
83
|
1
|
|
|
1
|
0
|
76
|
my $class = shift; |
84
|
1
|
|
|
|
|
11
|
my @pars = |
85
|
|
|
|
|
|
|
( |
86
|
|
|
|
|
|
|
'timecol', 1, 't' |
87
|
|
|
|
|
|
|
, 'Column for the (time) coordinate the ODE shall be advanced on.' |
88
|
|
|
|
|
|
|
. ' In the ODE, you can access it via [1] if the column is 1 (just like' |
89
|
|
|
|
|
|
|
. ' any other variable of the (interpolated) input data).' |
90
|
|
|
|
|
|
|
, 'ode', 'D0 = 1', 'e' |
91
|
|
|
|
|
|
|
, 'The ordinary differential equation system. The return value of the generated' |
92
|
|
|
|
|
|
|
. ' function does not matter, only that you set the values of the D array.' |
93
|
|
|
|
|
|
|
, 'varinit', [], 'i' |
94
|
|
|
|
|
|
|
, 'array of initial values; must match number of derivatives from ODE' |
95
|
|
|
|
|
|
|
, 'vartitle', [], '' |
96
|
|
|
|
|
|
|
, 'array of column titles for the integrated variables' |
97
|
|
|
|
|
|
|
, 'const', [], 'n' |
98
|
|
|
|
|
|
|
, 'array of constants to use in ODE' |
99
|
|
|
|
|
|
|
, 'rksteps', '4', 'k' |
100
|
|
|
|
|
|
|
, 'steps (stages) of the RK integration scheme, only 4 supported right now' |
101
|
|
|
|
|
|
|
, 'timestep', 0, 's' |
102
|
|
|
|
|
|
|
, 'desired time step size (see timediv)' |
103
|
|
|
|
|
|
|
, 'timediv', 10, '' |
104
|
|
|
|
|
|
|
, 'Divide input time intervals by that to get the integration time step. If' |
105
|
|
|
|
|
|
|
. ' timestep is set to non-zero, still an integer division is used, but one that' |
106
|
|
|
|
|
|
|
. ' yields a step close to the desired one (subject to rounding).' |
107
|
|
|
|
|
|
|
, 'interpolate', 'spline', '' |
108
|
|
|
|
|
|
|
, 'type of interpolation to use for intermediate points: spline or linear' |
109
|
|
|
|
|
|
|
, 'debug', 0, 'd' |
110
|
|
|
|
|
|
|
, 'give some info that may help debugging, >1 increasing verbosity' |
111
|
|
|
|
|
|
|
, 'plainperl', 0, '' |
112
|
|
|
|
|
|
|
, 'Use plain Perl syntax for formula for full force without confusing the' |
113
|
|
|
|
|
|
|
. ' intermediate parser.' |
114
|
|
|
|
|
|
|
); |
115
|
1
|
|
|
|
|
14
|
return $class->SUPER::new |
116
|
|
|
|
|
|
|
({ |
117
|
|
|
|
|
|
|
parconf=>{ info=>$infostring } |
118
|
|
|
|
|
|
|
, pardef=>\@pars |
119
|
|
|
|
|
|
|
, filemode=>1 |
120
|
|
|
|
|
|
|
, pipemode=>1 |
121
|
|
|
|
|
|
|
, pipe_init=>\&prepare |
122
|
|
|
|
|
|
|
, pipe_file=>\&process_file |
123
|
|
|
|
|
|
|
}); |
124
|
|
|
|
|
|
|
} |
125
|
|
|
|
|
|
|
|
126
|
|
|
|
|
|
|
sub prepare |
127
|
|
|
|
|
|
|
{ |
128
|
1
|
|
|
1
|
0
|
3
|
my $self = shift; |
129
|
1
|
|
|
|
|
3
|
my $p = $self->{param}; |
130
|
|
|
|
|
|
|
|
131
|
1
|
|
|
|
|
3
|
$p->{ode} = shift(@{$self->{argv}}) |
132
|
1
|
50
|
|
|
|
2
|
if(@{$self->{argv}}); |
|
1
|
|
|
|
|
3
|
|
133
|
|
|
|
|
|
|
$p->{varinit} = $self->{argv} |
134
|
1
|
50
|
|
|
|
3
|
if(@{$self->{argv}}); |
|
1
|
|
|
|
|
3
|
|
135
|
|
|
|
|
|
|
|
136
|
1
|
|
|
|
|
2
|
$self->{rk} = {}; |
137
|
1
|
|
|
|
|
3
|
my $rk = $self->{rk}; |
138
|
1
|
|
|
|
|
2
|
$rk->{stages} = 0; |
139
|
1
|
|
|
|
|
3
|
$rk->{a} = []; |
140
|
1
|
|
|
|
|
3
|
$rk->{b} = []; |
141
|
1
|
|
|
|
|
2
|
$rk->{c} = []; |
142
|
|
|
|
|
|
|
|
143
|
|
|
|
|
|
|
# The generic rksteps might be something funny (like 33, 45) in case |
144
|
|
|
|
|
|
|
# I intend to introduce methods that differ in stages and order. |
145
|
1
|
50
|
|
|
|
3
|
if($p->{rksteps} == 1) # Euler |
146
|
|
|
|
|
|
|
{ |
147
|
0
|
|
|
|
|
0
|
$rk->{stages} = 1; |
148
|
0
|
|
|
|
|
0
|
$rk->{a} = [ [0] ]; |
149
|
0
|
|
|
|
|
0
|
$rk->{b} = [ 1 ]; |
150
|
0
|
|
|
|
|
0
|
$rk->{c} = [ 0 ]; |
151
|
|
|
|
|
|
|
} |
152
|
1
|
50
|
|
|
|
3
|
if($p->{rksteps} == 2) # Heun |
153
|
|
|
|
|
|
|
{ |
154
|
0
|
|
|
|
|
0
|
$rk->{stages} = 2; |
155
|
|
|
|
|
|
|
$rk->{a} = |
156
|
|
|
|
|
|
|
[ |
157
|
0
|
|
|
|
|
0
|
[ 0, 0 ] |
158
|
|
|
|
|
|
|
, [ 1, 0 ] |
159
|
|
|
|
|
|
|
]; |
160
|
0
|
|
|
|
|
0
|
$rk->{b} = [ 0.5, 0.5 ]; |
161
|
0
|
|
|
|
|
0
|
$rk->{c} = [ 0, 1 ]; |
162
|
|
|
|
|
|
|
} |
163
|
1
|
50
|
|
|
|
4
|
if($p->{rksteps} == 3) # Simpson |
164
|
|
|
|
|
|
|
{ |
165
|
0
|
|
|
|
|
0
|
$rk->{stages} = 3; |
166
|
|
|
|
|
|
|
$rk->{a} = |
167
|
|
|
|
|
|
|
[ |
168
|
0
|
|
|
|
|
0
|
[ 0, 0, 0 ] |
169
|
|
|
|
|
|
|
, [ 0.5, 0, 0 ] |
170
|
|
|
|
|
|
|
, [ -1, 2, 0 ] |
171
|
|
|
|
|
|
|
]; |
172
|
0
|
|
|
|
|
0
|
$rk->{b} = [ 1/6, 4/6, 1/6 ]; |
173
|
0
|
|
|
|
|
0
|
$rk->{c} = [ 0, 0.5, 1 ]; |
174
|
|
|
|
|
|
|
} |
175
|
1
|
50
|
|
|
|
3
|
if($p->{rksteps} == 4) # RK44 method |
176
|
|
|
|
|
|
|
{ |
177
|
1
|
|
|
|
|
2
|
$rk->{stages} = 4; |
178
|
|
|
|
|
|
|
$rk->{a} = |
179
|
|
|
|
|
|
|
[ |
180
|
1
|
|
|
|
|
4
|
[0, 0, 0, 0] |
181
|
|
|
|
|
|
|
, [0.5, 0, 0, 0] |
182
|
|
|
|
|
|
|
, [0, 0.5, 0, 0] |
183
|
|
|
|
|
|
|
, [0, 0, 1, 0] |
184
|
|
|
|
|
|
|
]; |
185
|
1
|
|
|
|
|
3
|
$rk->{b} = [1./6, 1./3, 1./3, 1./6]; |
186
|
1
|
|
|
|
|
2
|
$rk->{c} = [0, 0.5, 0.5, 1]; |
187
|
|
|
|
|
|
|
} |
188
|
|
|
|
|
|
|
|
189
|
|
|
|
|
|
|
return $self->error("Invalid RK setup (nothing for $p->{rksteps} stages).") |
190
|
1
|
50
|
|
|
|
3
|
unless($rk->{stages}); |
191
|
1
|
50
|
|
|
|
10
|
if($p->{debug}) |
192
|
|
|
|
|
|
|
{ |
193
|
0
|
|
|
|
|
0
|
print STDERR "Using RK scheme with $rk->{stages} stages, tableau:\n"; |
194
|
0
|
|
|
|
|
0
|
for(my $s=0; $s<$rk->{stages}; ++$s) |
195
|
|
|
|
|
|
|
{ |
196
|
|
|
|
|
|
|
print STDERR sprintf( '%5.3f |'.(' %5.3f' x $rk->{stages})."\n" |
197
|
0
|
|
|
|
|
0
|
, $rk->{c}[$s], @{$rk->{a}[$s]} ); |
|
0
|
|
|
|
|
0
|
|
198
|
|
|
|
|
|
|
} |
199
|
0
|
|
|
|
|
0
|
print STDERR '------|'.('------' x $rk->{stages})."\n"; |
200
|
|
|
|
|
|
|
print STDERR sprintf( ' |'.(' %5.3f' x $rk->{stages})."\n" |
201
|
0
|
|
|
|
|
0
|
, @{$rk->{b}} ); |
|
0
|
|
|
|
|
0
|
|
202
|
|
|
|
|
|
|
} |
203
|
|
|
|
|
|
|
|
204
|
|
|
|
|
|
|
# The ODE stored as sub reference. Work arrays are V and D in addition |
205
|
|
|
|
|
|
|
# to A and C. |
206
|
|
|
|
|
|
|
$self->{ode} = formula_function( $p->{ode} |
207
|
|
|
|
|
|
|
, {plainperl=>$p->{plainperl}, verbose=>$p->{debug}} |
208
|
1
|
|
|
|
|
6
|
, 'V', 'D' ); |
209
|
|
|
|
|
|
|
return $self->error("Failed to compile your ODE.") |
210
|
1
|
50
|
|
|
|
8
|
unless defined $self->{ode}; |
211
|
|
|
|
|
|
|
|
212
|
1
|
|
|
|
|
5
|
return 0; |
213
|
|
|
|
|
|
|
} |
214
|
|
|
|
|
|
|
|
215
|
|
|
|
|
|
|
sub process_file |
216
|
|
|
|
|
|
|
{ |
217
|
1
|
|
|
1
|
0
|
2
|
my $self = shift; |
218
|
1
|
|
|
|
|
3
|
my $p = $self->{param}; |
219
|
1
|
|
|
|
|
3
|
my $txd = $self->{txd}; |
220
|
|
|
|
|
|
|
|
221
|
1
|
|
|
|
|
2
|
$self->{A} = []; |
222
|
1
|
|
|
|
|
3
|
$self->{C} = []; |
223
|
1
|
|
|
|
|
2
|
@{$self->{C}} = @{$p->{const}}; |
|
1
|
|
|
|
|
3
|
|
|
1
|
|
|
|
|
2
|
|
224
|
|
|
|
|
|
|
|
225
|
1
|
|
|
|
|
4
|
my $cols = $txd->columns(); |
226
|
1
|
50
|
33
|
|
|
4
|
unless($cols > 0 and @{$txd->{data}}) |
|
1
|
|
|
|
|
5
|
|
227
|
|
|
|
|
|
|
{ |
228
|
0
|
|
|
|
|
0
|
print STDERR "No data?\n"; |
229
|
0
|
|
|
|
|
0
|
$txd->write_all($self->{out}); |
230
|
0
|
|
|
|
|
0
|
return; |
231
|
|
|
|
|
|
|
} |
232
|
1
|
|
|
|
|
3
|
my $tc = $p->{timecol}-1; |
233
|
1
|
50
|
33
|
|
|
6
|
if($tc < 0 or $tc >= $cols) |
234
|
|
|
|
|
|
|
{ |
235
|
0
|
|
|
|
|
0
|
$txd->{data} = []; |
236
|
0
|
|
|
|
|
0
|
$txd->{titles} = []; |
237
|
0
|
|
|
|
|
0
|
print STDERR "Bad time index.\n"; |
238
|
0
|
|
|
|
|
0
|
$txd->write_all($self->{out}); |
239
|
0
|
|
|
|
|
0
|
return; |
240
|
|
|
|
|
|
|
} |
241
|
|
|
|
|
|
|
|
242
|
|
|
|
|
|
|
# The initial values tell us how many variables to expect, |
243
|
|
|
|
|
|
|
# prepare titles for added columns and also set initial |
244
|
|
|
|
|
|
|
# values. |
245
|
1
|
|
|
|
|
2
|
my $vars = @{$p->{varinit}}; |
|
1
|
|
|
|
|
3
|
|
246
|
1
|
|
|
|
|
2
|
my $vi=0; |
247
|
1
|
50
|
|
|
|
2
|
if(@{$txd->{raw_header}}){ $txd->write_header($self->{out}); } |
|
1
|
|
|
|
|
9
|
|
|
0
|
|
|
|
|
0
|
|
248
|
1
|
50
|
|
|
|
4
|
if(@{$txd->{titles}}) |
|
1
|
|
|
|
|
6
|
|
249
|
|
|
|
|
|
|
{ |
250
|
0
|
|
|
|
|
0
|
for(my $vi=0; $vi<$vars; ++$vi) |
251
|
|
|
|
|
|
|
{ |
252
|
|
|
|
|
|
|
$txd->{titles}[$cols+$vi] = defined $p->{vartitle}[$vi] |
253
|
0
|
0
|
|
|
|
0
|
? $p->{vartitle}[$vi] |
254
|
|
|
|
|
|
|
: 'var'.($vi+1); |
255
|
|
|
|
|
|
|
} |
256
|
0
|
|
|
|
|
0
|
print {$self->{out}} ${$txd->title_line()}; |
|
0
|
|
|
|
|
0
|
|
|
0
|
|
|
|
|
0
|
|
257
|
|
|
|
|
|
|
} |
258
|
|
|
|
|
|
|
# We'll print data on the fly, not bothering to clog memory with |
259
|
|
|
|
|
|
|
# storage of the integrated variables. Also might interfere with the |
260
|
|
|
|
|
|
|
# interpolation otherwise. |
261
|
1
|
|
|
|
|
3
|
my @val = @{$p->{varinit}}; |
|
1
|
|
|
|
|
5
|
|
262
|
1
|
|
|
|
|
4
|
print {$self->{out}} ${$txd->data_line([ |
|
1
|
|
|
|
|
4
|
|
263
|
1
|
|
|
|
|
5
|
@{$txd->{data}[0]}[0..($cols-1)] |
|
1
|
|
|
|
|
15
|
|
264
|
|
|
|
|
|
|
, @val ])}; |
265
|
1
|
|
|
|
|
5
|
for(my $mi=1; $mi< @{$txd->{data}}; ++$mi) |
|
6
|
|
|
|
|
29
|
|
266
|
|
|
|
|
|
|
{ |
267
|
|
|
|
|
|
|
# Integrate from to to t1, using a fixed step that fits into the interval. |
268
|
5
|
|
|
|
|
17
|
my $t0 = $txd->{data}[$mi-1][$tc]; |
269
|
5
|
|
|
|
|
11
|
my $t1 = $txd->{data}[$mi][$tc]; |
270
|
|
|
|
|
|
|
my $div = $p->{timestep} |
271
|
|
|
|
|
|
|
? int(abs(($t1-$t0)/$p->{timestep})+0.5) |
272
|
5
|
50
|
|
|
|
19
|
: $p->{timediv}; |
273
|
5
|
50
|
|
|
|
17
|
$div = 1 |
274
|
|
|
|
|
|
|
if $div < 1; |
275
|
5
|
|
|
|
|
14
|
my $step = ($t1-$t0)/$div; |
276
|
|
|
|
|
|
|
print STDERR "int $t0 to $t1 div $div step $step\n" |
277
|
5
|
50
|
|
|
|
15
|
if $p->{debug}; |
278
|
5
|
|
|
|
|
14
|
for(my $si=0; $si<$div; ++$si) |
279
|
|
|
|
|
|
|
{ |
280
|
50
|
|
|
|
|
154
|
$self->rk_step($t0+$si*$step, $step, \@val); |
281
|
|
|
|
|
|
|
} |
282
|
5
|
|
|
|
|
12
|
print {$self->{out}} ${$txd->data_line([ |
|
5
|
|
|
|
|
16
|
|
283
|
5
|
|
|
|
|
15
|
@{$txd->{data}[$mi]}[0..($cols-1)] |
|
5
|
|
|
|
|
53
|
|
284
|
|
|
|
|
|
|
, @val ])}; |
285
|
|
|
|
|
|
|
} |
286
|
|
|
|
|
|
|
} |
287
|
|
|
|
|
|
|
|
288
|
|
|
|
|
|
|
# Compute one RK step with given time increment. |
289
|
|
|
|
|
|
|
sub rk_step |
290
|
|
|
|
|
|
|
{ |
291
|
50
|
|
|
50
|
0
|
85
|
my $self = shift; |
292
|
50
|
|
|
|
|
99
|
my ($t, $dt, $val) = @_; |
293
|
|
|
|
|
|
|
|
294
|
50
|
|
|
|
|
89
|
my $rk = $self->{rk}; |
295
|
50
|
|
|
|
|
72
|
my $vars = @{$val}; |
|
50
|
|
|
|
|
85
|
|
296
|
|
|
|
|
|
|
|
297
|
50
|
|
|
|
|
91
|
my @work; # Storage for the derivatives in the stages. |
298
|
|
|
|
|
|
|
my @tmp; # Storage for the current variables. |
299
|
|
|
|
|
|
|
|
300
|
|
|
|
|
|
|
# Initialise them with the correct size. |
301
|
50
|
|
|
|
|
120
|
for(my $s=0; $s<$rk->{stages}; ++$s) |
302
|
|
|
|
|
|
|
{ |
303
|
200
|
|
|
|
|
376
|
for(my $v=0; $v<$vars; ++$v) |
304
|
|
|
|
|
|
|
{ |
305
|
400
|
|
|
|
|
962
|
$work[$s][$v] = 0; |
306
|
|
|
|
|
|
|
} |
307
|
|
|
|
|
|
|
} |
308
|
50
|
|
|
|
|
109
|
for(my $v=0; $v<$vars; ++$v) |
309
|
|
|
|
|
|
|
{ |
310
|
100
|
|
|
|
|
196
|
$tmp[$v] = 0; |
311
|
|
|
|
|
|
|
} |
312
|
|
|
|
|
|
|
|
313
|
|
|
|
|
|
|
# Collect the stage derivatives. |
314
|
50
|
|
|
|
|
137
|
$self->eval_ode($t, $val, $work[0]); |
315
|
0
|
|
|
|
|
0
|
print STDERR "deriv 0: @{$work[0]}\n" |
316
|
50
|
50
|
|
|
|
147
|
if $self->{param}{debug} > 1; |
317
|
50
|
|
|
|
|
129
|
for(my $stage=1; $stage < $rk->{stages}; ++$stage) |
318
|
|
|
|
|
|
|
{ |
319
|
150
|
|
|
|
|
295
|
for(@tmp){ $_ = 0 } |
|
300
|
|
|
|
|
517
|
|
320
|
150
|
|
|
|
|
337
|
for(my $substage = 0; $substage < $stage; ++$substage) |
321
|
|
|
|
|
|
|
{ |
322
|
300
|
100
|
|
|
|
720
|
if($rk->{a}[$stage][$substage] != 0) |
323
|
|
|
|
|
|
|
{ # Does the condition really save work? |
324
|
150
|
|
|
|
|
294
|
for(my $i=0; $i<$vars; ++$i) |
325
|
|
|
|
|
|
|
{ |
326
|
300
|
|
|
|
|
898
|
$tmp[$i] += $rk->{a}[$stage][$substage]*$work[$substage][$i]; |
327
|
|
|
|
|
|
|
} |
328
|
|
|
|
|
|
|
} |
329
|
|
|
|
|
|
|
} |
330
|
150
|
|
|
|
|
333
|
for(my $i=0; $i<$vars; ++$i){ $tmp[$i] *= $dt; $tmp[$i] += $val->[$i]; } |
|
300
|
|
|
|
|
446
|
|
|
300
|
|
|
|
|
667
|
|
331
|
150
|
|
|
|
|
534
|
$self->eval_ode($t+$rk->{c}[$stage]*$dt, \@tmp, $work[$stage]); |
332
|
0
|
|
|
|
|
0
|
print STDERR "deriv $stage: @{$work[$stage]}\n" |
333
|
150
|
50
|
|
|
|
624
|
if $self->{param}{debug} > 1; |
334
|
|
|
|
|
|
|
} |
335
|
|
|
|
|
|
|
|
336
|
|
|
|
|
|
|
# Compute the definite derivative, apply and be done. |
337
|
50
|
|
|
|
|
102
|
for(@tmp){ $_ = 0 } |
|
100
|
|
|
|
|
180
|
|
338
|
50
|
|
|
|
|
117
|
for(my $stage=0; $stage < $rk->{stages}; ++$stage) |
339
|
|
|
|
|
|
|
{ |
340
|
200
|
|
|
|
|
394
|
for(my $i=0; $i<$vars; ++$i) |
341
|
|
|
|
|
|
|
{ |
342
|
400
|
|
|
|
|
979
|
$tmp[$i] += $rk->{b}[$stage]*$work[$stage][$i]; |
343
|
|
|
|
|
|
|
} |
344
|
|
|
|
|
|
|
} |
345
|
50
|
|
|
|
|
117
|
for(my $i=0; $i<$vars; ++$i) |
346
|
|
|
|
|
|
|
{ |
347
|
100
|
|
|
|
|
377
|
$val->[$i] += $tmp[$i]*$dt; |
348
|
|
|
|
|
|
|
} |
349
|
|
|
|
|
|
|
} |
350
|
|
|
|
|
|
|
|
351
|
|
|
|
|
|
|
# Evaluate the ODE once, interpolating in the input data for time-varying |
352
|
|
|
|
|
|
|
# parameters. |
353
|
|
|
|
|
|
|
sub eval_ode |
354
|
|
|
|
|
|
|
{ |
355
|
200
|
|
|
200
|
0
|
346
|
my $self = shift; |
356
|
200
|
|
|
|
|
372
|
my ($t, $var, $deriv) = @_; |
357
|
200
|
|
|
|
|
284
|
my @fd; # interpolated data set here |
358
|
200
|
|
|
|
|
603
|
$fd[0] = $self->{txd}->set_of($t, $self->{param}{timecol}-1); |
359
|
0
|
|
|
|
|
0
|
print STDERR "fd: @{$fd[0]}\n" |
360
|
200
|
50
|
|
|
|
555
|
if $self->{param}{debug} > 1; |
361
|
0
|
|
|
|
|
0
|
print STDERR "V: @{$var}\n" |
362
|
200
|
50
|
|
|
|
426
|
if $self->{param}{debug} > 1; |
363
|
|
|
|
|
|
|
# @{$deriv} == 0 since rk_step intialised it |
364
|
200
|
|
|
|
|
5079
|
$self->{ode}->(\@fd, $self->{A}, $self->{C}, $var, $deriv); |
365
|
|
|
|
|
|
|
} |
366
|
|
|
|
|
|
|
|
367
|
|
|
|
|
|
|
1; |