Math::ODE - Solve N-th Order Ordinary Differential Equations

This module allows you to solve N-th Order Ordinary Differential Equations with as little pain as possible. Currently, only IVP's (initial value problems) are supported, but native support for BVP's (boundary value problems) may be added in the future. To solve N-th order equations, you must first turn it into a system of N first order equations, as in MATLAB.

	use Math::ODE;
	# create new object that stores data in a file 
	# and solve the given equation(s) on the interval
	# [0,10], with initial condition y(t0) = 0
	my $o = new Math::ODE ( file => '/home/user/data',
                        step => 0.1,
                        initial => [0], 
                        ODE => [ \&DE1 ], 
                        t0 => 1,
                        tf => 10 );
	$o->evolve();
	# solve the equation y' = 1/$t
	# $t is the independent variable, a scalar
	# $y is the dependent variable, an array reference
	sub DE1 { my ($t,$y) = @_; return 1/$t; }

• $o->evolve()
  
  Evolves the equations from $o->t0 to $o->tf using a 4th Order Classic Runge-Kutta method.
  
  

  	# Example 1: Solve  y'' + y = 0, y(0) = 0, y'(0) = 1
  	#            Solution: y = sin(x)
  
  

  
  
  

  	use Math::ODE;
  	my $o = new Math::ODE ( file => 'data',
  				verbose => 1,
  				step => 0.1,
  				initial => [0,1], 
  				ODE => [ \&DE0, \&DE1 ], 
  				t0 => 0,
  				tf => 10 );
  	$o->evolve;
  
  

  
  
  

  	# to plot data in gnuplot: plot 'data' using 1:2, sin(x)
  
  

  
  
  

  	# y'' + y = 0
  	sub DE0 { my ($t,$y) = @_; return $y->[1]; }
  	sub DE1 { my ($t,$y) = @_; return -$y->[0]; }
  
  

  
  
  To turn y'' + y = 0 into a system, we will imagine a vector with two components, called y. Now let the first component y0 = y and the second component y1 = y0' . Now rewrite the equation in terms of these variables. It will be y1' + y0 = 0. Solving for y1' we get y1' = -y0. Now the vector y' has compenents y0' = y1 and y1' = -y0 . These are the equations we put in our DE* sub's, except we use an arrayref. The DE0 sub corresponds to our y0' component and the DE1 sub corresponds to our y1' component. This can be *ahem* easily generalized to any N-th Oder equation or system of equations. An example of a system of second order equations is given in example/de4. Please look in there for other examples as well.
  
  The initial arrayref corresponds to the initial conditions of each of the components of the dependent variable vector. [0,1] means that the first component at $o->t0 (which happens to be 0) will have the value 0 and the second component at $o->t0 (which is 10) will have the value 1.
• $o->ODE($F)
  
  Sets the equations to be solved to $F. $F must be an arrayref of coderefs and the same length as $F, or Bad Things May Occur (tm).
• $o->initial($arrayref)
  
  Set the initial conditions to $arrayref. Returns an arrayref of initial conditions if no arguments are given. Must be the same size as $o->ODE or autovivisection will ensue.
• $o->step($s)
  
  This will set the step size (length of the segments of the interval) to $s, which must be positive. The default is 0.01. Returns the step size if no arguments are given.
• $o->t0($a)
  
  This will set the initial (leftmost) point on the interval to $a. Returns the initial point on the interval if no arguments are given.
• $o->tf($b)
  
  This will set the last (rightmost) point on the interval to $b. Returns the end point on the interval if no arguments are given.
• $o->file($somefile)
  
  Save data in $somefile. Returns the file in which the data is being saved if no arguments are given. If no file is specified, data is printed to STDOUT. The data file will have $N+1 columns (where $N is the number of equations to solve), and can be fed directly to gnuplot. The first column is the independent variable, and the remaining are the first through nth components of the dependent vector. Examples of graphing the data file are in the example/ directory of the source distribution.
• $o->verbose($number)
  
  Sets the verbosity of debugging output. The default of 1 will currently only cause evolve() to warn you when it gets NaN or Inf, instead of silently returning undef, when verbosity is set to 0. Setting the verbosity to 0 is useful if you are trying to write code to solve boundary value problems with the shooting method. You will be guessing initial conditions, and don't feel like getting a warning on every guess that blows up (which will be many.) Setting the verbosity to 2 will cause a message like the following:
  
  

  	Exiting RK4 with t=9.9 ,$y1 = ['-0.544013766248772833','-0.839075464413064726'];
  
  

  
  
  to be printed on every increment of the independent variable $t. These are the values that the 4th Order Runge-Kutta returned for the current value of $t.
• >Returns the maximum error from the computed values and a reference to a list of code references.
  
  

Jonathan Leto <jonathan@leto.net>

Boyce, DiPrima "Elementary Differential Equations" 5th Ed.

Orwant, Hietaniemi, Macdonald "Mastering Algorithms with Perl" Ch. 16.

Copyright (c) 2001-2008 by Jonathan Leto. All rights reserved.

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