File Coverage

blib/lib/Math/ODE.pm

Criterion	Covered	Total	%
statement	101	112	90.1
branch	22	36	61.1
condition	2	7	28.5
subroutine	12	14	85.7
pod	0	7	0.0
total	137	176	77.8

package Math::ODE;
use strict;

use Data::Dumper;
use Carp;
my $VERSION = '0.05_01';

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

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

    # don't clobber the initial condition in case we want to do multiple runs
    # with the same object
    @$y = @$initial;

    if( defined $file ){
            open($fh,'>', $file) or croak "$file: $!";
    }

    $self->_clear_values;

    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;

        my $str = join $delim,  map { sprintf "%0.12f", $_ } ($t, @$y);
        chop $str;

        if( defined $file ){
            print $fh "$str\n";
        } elsif ( ! $self->{keep_values} ) {
            print "$str\n";
        }
    }
    close $fh if defined $file;
    return $self;
}

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,$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);
        my @indices = ( 0 .. $self->{N} - 1 );

        map { $w1[$_]  = $h * &{ $F->[$_] }($t,$y)           } @indices;
        map { $q->[$_] = $y->[$_] + 0.5*$w1[$_]              } @indices;

        map { $w2[$_]  = $h * &{ $F->[$_] }($t + 0.5*$h,$q)  } @indices;
        map { $q->[$_] = $y->[$_] + 0.5*$w2[$_]              } @indices;

        map { $w3[$_]  = $h * &{ $F->[$_] }($t + 0.5*$h,$q)  } @indices;
        map { $q->[$_] = $y->[$_] + $w3[$_]                  } @indices;

        map { $w4[$_]  = $h * &{ $F->[$_] }($t + $h,$q)      } @indices;
        map { $y->[$_]+= ( $w1[$_] + 2 * $w2[$_] + 2 * $w3[$_] + $w4[$_])/6 } @indices;

        $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 step
{
    my ($self,$step) = @_;
    if (defined $step){
        croak "Stepsize must be strictly between zero and one"
                if ($step <= 0 || $step >= 1);

        $self->{step} = $step;
    }
    return $self->{step};
}

sub initial
{
    my ($self, $initial) = @_;

    croak "not an array ref" unless ref $initial eq "ARRAY";

    $self->{initial} = $initial;

}

sub _clear_values
{
    my $self = shift;
    $self->{values} = {}
}

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

    if( $self->{N} != scalar(  @{ $args{initial } }) ){
                croak "Must have same number of initial conditions as equations!";
    }

    croak "Stepsize must be positive!"               if $self->{step} <= 0;
    croak "\$self->t0 must be less than \$self->tf!" if $self->{t0}   >= $self->{tf};

    return $self;
}

sub max_error
{
    my ($self, $sols) = @_;

    my $max_error = 0;

    for my $pt ( sort keys %{$self->{values}} ){
        my $k = 0;
        my @vals = $self->values_at($pt);

        for my $sol ( @$sols ) {
            my $res  = abs( $vals[$k]  - &$sol($pt) );
            $max_error = $res if ($res > $max_error);
            $k++;
        }
    }
    $max_error;
}

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

42;
__END__