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

Math-ODE documentation Contained in the Math-ODE distribution.

Index

Code Index:

NAME

Top

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

DESCRIPTION

Top

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.

SYNOPSIS

Top

	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; }

AUTHOR

Top

Jonathan Leto <jonathan@leto.net>

SEE ALSO

Top

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

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

COPYRIGHT

Top

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

LICENSE AGREEMENT

Top

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

Math-ODE documentation Contained in the Math-ODE distribution. 
package Math::ODE;
use strict;

require 5.005;
use Data::Dumper;
use Carp;
my $VERSION = '0.04';

$Data::Dumper::Varname = "y";
$Data::Dumper::Indent = 0;

sub evolve {
	my $self = shift;
	my ($F,$h,$t,$y,$file) = map{ $self->{$_} } qw(ODE step t0 initial file);
	my $delim = $self->{csv} ? ',' : ($self->{delim} || $self->{delimeter} || " ");
	
	if( defined $file ){
        	open(FD, ">$file") or croak "$file: $!";
	}

    while ( $t < $self->{tf} ){
        # use Runge Kutta 4th order to step from $t to $t + $h
        $y = _RK4($self,$t,$y);

        warn "Exiting RK4 with t=$t ," . Dumper($y) . "\n" if( $self->{verbose} > 1 );

        for my $i ( 0 .. $self->{N}-1 ){
            # check for under/over flow
            next unless $y->[$i] =~ qr/nan|infinity/i;
            warn "Bailing out, over/under flow at t=$t,y->[$i] = $y->[$i]" if $self->{verbose};
            return undef;
        }
        $t += $h;

        if( defined $file ){
            my $str = join $delim,  map { sprintf "%0.12f", $_ } ($t, @$y);
		    chop $str;
            print FD "$str\n";
	    }
    }
	close FD if defined $file;
    return 42;
}

sub _RK4 {
        # $t = dependent variable
        # $y = $N - vector of independent variables
        # $h = step size
        # $F = arrayref of coderefs of the equations to solve
        my $self = shift;
        my ($t, $y) = @_;
        my $F = $self->{ODE};
        my $h = $self->{step};

        ## w vectors hold constants for equations
        ## each $q holds a modified $y vector to feed to the next
        ## for loop ( ie $y + $w1/2 , etc ... )
        my (@w1,@w2,@w3,@w4,$q,$i);

        for $i ( 0 .. $self->{N}-1 ){ $w1[$i]  = $h * &{ $F->[$i] }($t,$y);           }
        for $i ( 0 .. $self->{N}-1 ){ $q->[$i] = $y->[$i] + 0.5*$w1[$i];              }

        for $i ( 0 .. $self->{N}-1 ){ $w2[$i]  = $h * &{ $F->[$i] }($t + 0.5*$h,$q);  }
        for $i ( 0 .. $self->{N}-1 ){ $q->[$i] = $y->[$i] + 0.5*$w2[$i];              }

        for $i ( 0 .. $self->{N}-1 ){ $w3[$i]  = $h * &{ $F->[$i] }($t + 0.5*$h,$q);  }
        for $i ( 0 .. $self->{N}-1 ){ $q->[$i] = $y->[$i] + $w3[$i];                  }

        for $i ( 0 .. $self->{N}-1 ){ $w4[$i]  = $h * &{ $F->[$i] }($t + $h,$q);      }


        for $i ( 0 .. $self->{N}-1 ){ $y->[$i] += ( $w1[$i] + 2 * $w2[$i] + 2 * $w3[$i] + $w4[$i])/6; }

	    $self->_store_values( $t + $h, $y );
	
        return $y;
}
sub _store_values {
	my ($self,$t, $y) = @_;
	return unless  $self->{keep_values};
	my $s = sprintf '%0.12f', $t ; 
	push @{ $self->{values}{$s} }, @$y;
}
sub values_at {
	my ($self,$t, %args) = @_;
    if ($self->{keep_values}){
	    return @{ $self->{values}{sprintf('%0.12f',$t)} };
    } else {
        warn "Values were not kept because keep_values was set to 0";
        return;
    }
}
# because Math::ODE implements a 4th order Runge-Kutta method
sub error {  $_[0]->{step} ** 4 }

sub _init {
	my ($self,%args) = @_;

	# defaults
	$self->{keep_values} = 1;
	$self->{verbose}     = 1;
	$self->{step}        = 0.1; 
	$self->{csv}         = 0;
	$self->{N}           = scalar( @{ $args{ODE} } ) || 1;

	@$self{keys %args} = values %args;
	$self->{values} = {};

    if( $self->{N} != scalar(  @{ $args{initial } }) ){
                croak "Must have same number of initial conditions as equations!";
    }
	if( $self->{step} <= 0  ){
		croak "Stepsize must be positive!";
	}
	if( $self->{t0} >= $self->{tf} ){
		croak "\$self->t0 must be less than \$self->tf!";
	}
    return $self;
}
sub max_error {
    my ($self, $sols) = @_;

    my $max_error = 0;

    for my $pt ( sort keys %{$self->{values}} ){
        my $k = 0;
        for my $sol ( @$sols ) {
            my @vals = $self->values_at($pt);
            my $res  = abs( $vals[$k]  - &$sol($pt) );
            $max_error = $res if ($res > $max_error);
            print "pt=$pt, res=$res\n" if ($res > $self->error && debug() );
            $k++;
        }
    }
    $max_error;
}
sub debug { 0 }

sub new {
	my $class = shift;
	my $self = {};
	bless($self, $class);
	$self->_init(@_);
}
42;
__END__