Hello all,
Just in case it's of interest to anyone else, I've implemented a
"cumuintover" function. Analogously to "cumusumover" and
"cumuprodover" it returns a piddle with the same dimensions as the
input containing the integral over the first dimension.
Note that it will crash if you feed it a piddle with a first dimension
of less than 6. If it's of use to anyone, I could modify it to behave
like intover and work for dims < 6.
This is my first use of PP, so I welcome any feedback folks might
have.
Attached is a demo of the code using inline, along with the pp_def.
Here's what I get:
perl demo_cumuintover.pl
x | num | ana | rel err
0 0 0 nan
1 2.13 2.13 8.43e-08
2 14.1 14.1 3.68e-08
3 56.5 56.5 1.91e-08
4 162 162 1.23e-08
5 377 377 9.35e-09
6 758 758 7.95e-09
7 1.38e+03 1.38e+03 7.37e-09
8 2.33e+03 2.33e+03 7.3e-09
9 3.7e+03 3.7e+03 7.62e-09
10 5.63e+03 5.63e+03 8.3e-09
-roban
#!/usr/bin/perl -w
use strict;
use PDL;
use PDL::NiceSlice;
use Inline qw(Pdlpp);
my $dx = 0.1;
my $x = $dx * pdl [0..100];
my $y = $x**2 + 2*$x**3 + exp($x/2);
my $g_num1 = $dx * PDL::cumuintover($y);
my $g_ana = $x**3/3 + 2*$x**4/4 + 2 * exp($x/2) - 2;
print " x | num | ana | rel err\n";
wcols "%8.3g ", $x(0:-1:10), $g_num1(0:-1:10), $g_ana(0:-1:10),
($g_num1(0:-1:10) - $g_ana(0:-1:10))/$g_ana(0:-1:10);
no PDL::NiceSlice;
__DATA__
__Pdlpp__
# from equation 4.1.14, Press et al. 2nd Ed chapter 4.1
pp_def('cumuintover',
Pars => 'f(n); int+ [o]g(n);',
Code => '
int i;
int ns = $SIZE(n);
$GENERIC(g) tmp = 0;
threadloop %{
tmp = (3.0/8.0) * $f(n=>0) +
(7.0/6.0) * $f(n=>1) +
(23.0/24.0) * $f(n=>2);
for (i=5;i<ns;i++) { $g(n=>i) = tmp +
(23./24.) * $f(n=>i-2) +
(7./6.) * $f(n=>i-1) +
(3./8.) * $f(n=>i);
tmp += $f(n=>i-2);
}
$g(n=>0) = 0;
$g(n=>1) = ($f(n=>0) + $f(n=>1))/2.0;
$g(n=>2) = ($f(n=>0) + 4. * $f(n=>1) + $f(n=>2) ) / 3.0;
$g(n=>3) = (3/8) * ($f(n=>0) + $f(n=>3) + 3. * ( $f(n=>1) + $f(n=>2) ) );
$g(n=>4) = (1./45.) * ( 14. * ( $f(n=>0) + $f(n=>4) )
+ 64. * ( $f(n=>1) + $f(n=>3) )
+ 24. * $f(n=>2)
);
%}
'
);
_______________________________________________
Perldl mailing list
[email protected]
http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
