View on
Dragos Constantinescu > PDL-RungeKutta > PDL::RungeKutta



Annotate this POD


New  1
Open  0
View/Report Bugs
Module Version: 0.01   Source  


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                                            

  # integration                                                                
  ($evt,$evy,$evd,$i,$j) =                                                     

  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 ^


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

* $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.
* $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.

* $t initial moment
* $yi piddle containing the initial conditions
* $h step
* $DE reference to differential equation
* $y piddle containing the result
* $del piddle containing the errors


Dragos Constantinescu <>



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


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: