Dragos Constantinescu >
PDL-RungeKutta >
PDL::RungeKutta

Module Version: 0.01
PDL::RungeKutta - Solve N-th order M dimensional ordinary differential equations using adaptive stepsize Runge Kutta method.

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.

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.

($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

($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

Dragos Constantinescu <dragos@venus.nipne.ro>

Math::ODE(3)

"Numerical Recipes in Fortran 77: The Art of Scientific Computing" Ch. 16.

http://lib-www.lanl.gov/numerical/index.html

Copyright (c) 2003 by Dragos Constantinescu. All rights reserved.

This package is free software; you can redistribute it and/or modify it under the same terms as Perl itself.

syntax highlighting: