#! /usr/bin/perl ######################################################################### # This Perl script is Copyright (c) 2002, Peter J Billam # # c/o P J B Computing, www.pjb.com.au # # # # This script is free software; you can redistribute it and/or # # modify it under the same terms as Perl itself. # ######################################################################### use Math::RungeKutta qw(:ALL); use Test::Simple tests => 6; my $func_evals = 0; my $i_test = 0; my $n_failed = 0; my $n_passed = 0; sub dydt { my ($t, @y) = @_; my @dydt; $dydt[$[] = $y[$[+1]; $dydt[$[+1] = 0.0 - $y[$[]; $func_evals++; return @dydt; } $twopi = 2.0 * 3.141592653589; %passmark0 = ( rk2=>0.2, rk4=>.0004, rk4_classical=>.0015, rk4_ralston=>.0015, epsilon=>.0001, errors=>.0003, ); %passmark1 = ( rk2=>0.04, rk4=>.00004, rk4_classical=>.0006, rk4_ralston=>.0006, epsilon=>.00001, errors=>.00001, ); foreach $algorithm ('rk2','rk4','rk4_classical','rk4_ralston') { $i_test++; $n = 16; $dt= $twopi / $n; @y = (0,1); $t=0; foreach (1..$n) { ($t, @y) = &{$algorithm}( \@y, \&dydt, $t, $dt ); } my $err0 = abs $y[$[]; my $err1 = abs ($y[$[+1]-1.0); ok (($err0 < $passmark0{$algorithm} && $err1 < $passmark1{$algorithm}), $algorithm); } $algorithm = 'rk4_auto'; my ($t_midpoint, @y_midpoint); MODE: foreach $mode ('epsilon','errors') { $i_test++; @y = (0,1); $t=0; $dt = 0.1; my $i = 0; my $epsilon; if ($mode eq 'epsilon') { $epsilon = .0001; } else { @errors = (.01, .0001); $epsilon = \@errors; } $func_evals = 0; while ($t+$dt < $twopi) { $i++; ($t, $dt, @y) = &rk4_auto( \@y, \&dydt, $t, $dt, $epsilon ); ($t_midpoint, @y_midpoint) = &rk4_auto_midpoint(); if ($func_evals > 500) { ok(0, "rk4_auto($mode) failed, $func_evals func evals"); next MODE; } } $i++; $dt = $twopi-$t; ($t, @y) = &rk4( \@y, \&dydt, $t, $dt ); my $err0 = abs $y[$[]; my $err1 = abs ($y[$[+1]-1.0); ok (($err0 < $passmark0{$mode} && $err1 < $passmark1{$mode}), "rk4_auto($mode)"); } __END__ =pod =head1 NAME test.pl - Perl script to test Math::RungeKutta.pm =head1 SYNOPSIS perl test.pl =head1 DESCRIPTION This script tests Math::RungeKutta.pm =head1 AUTHOR Peter J Billam http://www.pjb.com.au/comp/contact.html =head1 SEE ALSO Math::RungeKutta.pm , http://www.pjb.com.au/ , perl(1). =cut