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.