File Coverage

blib/lib/Text/NumericData/App/txdodeint.pm
Criterion Covered Total %
statement 117 154 75.9
branch 22 44 50.0
condition 2 6 33.3
subroutine 8 8 100.0
pod 0 5 0.0
total 149 217 68.6


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;