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.