NAME
PDL::RungeKutta - Solve N-th order M dimensional ordinary differential equations using adaptive stepsize Runge Kutta method.
DESCRIPTION
This module allows to solve N-th order M dimensional ordinary differential equations. It uses the adaptive stepsize control for fifth order Cash-Karp Runge-Kutta algorithm described in "Numerical Recipes in Fortran 77: The Art of Scientific Computing" Ch. 16.2. The errors are estimated as the difference between fifth order results and the embeded forth order results. To solve N-th order equations, you must first turn it into a system of N first order equations.
SYNOPSIS
Example: Solve y'' + y = 0, y(0) = 0, y'(0) = 1
use PDL;
use PDL::Math;
use PDL::NiceSlice;
use PDL::RungeKutta;
# y'' + y = 0, Solution: y = sin(t)
$Y0 = pdl(0,1); # y(0)=0, y'(0)=1 ( Y0=(f,g), f=y, g=y' )
@esargs=(); # extra arguments for error eval function
$t0 = 0; # initial moment
$dt0 = 0.1; # initial time step
$t1 = 1 0; # final moment
$eps = 1.e-6; # error
$verbose=1;
# integration
($evt,$evy,$evd,$i,$j) =
rkevolve($t0,$Y,$dt0,\&DE,$t1,$eps,\&error,\@esargs,$verbose);
$check=sin($evt);
wcols $evt,$evy((0)),$check,'test.dat';
sub DE { # differential eq
my ($t,$y)= @_;
my $yd=zeroes(2); # Y' ( = (f',g') = (y',y'') )
$yd(0).=$y(1); # f'=g ( = y' )
$yd(1).=-$y(0); # g'=-f ( =-y )
return $yd;
}
sub error { # error scale
my ($t,$Y) = @_;
my $es=ones(2); # constant scale
return $es;
}
Please see the other examples also.
Exported Functions
rkevolve
($t,$Y,$d,$i,$j) = rkevolve($t0,$Y0,$dt0,$DE,$t1,$eps,$erfcn,$efarg,$verbose)
This will solve the differential equation for DE with the initial conditions Y0 between t0 and t1 using adaptive step size control for Runge-Kutta.
- Input
- * $t0 scalar, initial moment
- * $Y0 one dimensional piddle wich contains the initial conditions. Number of elements is NxM.
- * $dt0 scalar, initial time step
- * $DE reference to the function for the differential equation. Please see Math::ODE(3) for a more detailed description on how to construct the equations. It should return a piddle with the same dimensions as $Y0.
- * $t1 scalar, final moment
- * $eps scalar, the requested error. Upon scaling with the output of erfcn this will give the requested error for each element of $Y.
- * $erfcn reference to the error scaling function. This function should return a pidle containing scaling factors for the requested error for each element of $Y. Please see Num. Rec. for more details.
- * $efarg reference to an array containing supplementary arguments for erfcn.
- * $verbose scalar, integer. If set to 1 details about progress are printed during calculation.
- Output
- * $t piddle containing the independent variable
- * $Y piddle containing the results
- * $d piddle containing the errors as estimated by Cash-Karp Runge-Kutta algorithm
- * $i total number of iterations
- * $j number of resets made in order to decrease the error
rk5
($y,$del) = rk5($t,$yi,$h,$DE)
This will carry out one Cash-Karp Runge-Kutta step. Could be useful if you want to do your own step size control or you just want equal steps.
- Input
- * $t initial moment
- * $yi piddle containing the initial conditions
- * $h step
- * $DE reference to differential equation
- Output
- * $y piddle containing the result
- * $del piddle containing the errors
AUTHOR
Dragos Constantinescu <dragos@venus.nipne.ro>
SEE ALSO
Math::ODE(3)
"Numerical Recipes in Fortran 77: The Art of Scientific Computing" Ch. 16.
http://lib-www.lanl.gov/numerical/index.html
COPYRIGHT
Copyright (c) 2003 by Dragos Constantinescu. All rights reserved.
LICENSE AGREEMENT
This package is free software; you can redistribute it and/or modify it under the same terms as Perl itself.
24 POD Errors
The following errors were encountered while parsing the POD:
- Around line 145:
'=item' outside of any '=over'
- Around line 147:
Expected text after =item, not a bullet
- Around line 149:
Expected text after =item, not a bullet
- Around line 152:
Expected text after =item, not a bullet
- Around line 154:
Expected text after =item, not a bullet
- Around line 158:
Expected text after =item, not a bullet
- Around line 160:
Expected text after =item, not a bullet
- Around line 163:
Expected text after =item, not a bullet
- Around line 167:
Expected text after =item, not a bullet
- Around line 170:
Expected text after =item, not a bullet
- Around line 175:
Expected text after =item, not a bullet
- Around line 177:
Expected text after =item, not a bullet
- Around line 179:
Expected text after =item, not a bullet
- Around line 182:
Expected text after =item, not a bullet
- Around line 184:
Expected text after =item, not a bullet
- Around line 186:
You forgot a '=back' before '=head2'
- Around line 193:
'=item' outside of any '=over'
- Around line 195:
Expected text after =item, not a bullet
- Around line 197:
Expected text after =item, not a bullet
- Around line 199:
Expected text after =item, not a bullet
- Around line 201:
Expected text after =item, not a bullet
- Around line 205:
Expected text after =item, not a bullet
- Around line 207:
Expected text after =item, not a bullet
- Around line 209:
You forgot a '=back' before '=head1'