File Coverage

blib/lib/Math/RungeKutta.pm

Criterion	Covered	Total	%
statement	464	522	88.8
branch	75	100	75.0
condition	12	27	44.4
subroutine	9	12	75.0
pod	6	9	66.6
total	566	670	84.4

# Math::RungeKutta.pm
#########################################################################
#        This Perl module is Copyright (c) 2002, Peter J Billam         #
#               c/o P J B Computing, www.pjb.com.au                     #
#                                                                       #
#     This module is free software; you can redistribute it and/or      #
#            modify it under the same terms as Perl itself.             #
#########################################################################

package Math::RungeKutta;
no strict; no warnings;
$VERSION = '1.07';
# gives a -w warning, but I'm afraid $VERSION .= ''; would confuse CPAN
require Exporter;
@ISA = qw(Exporter);
@EXPORT = qw(rk2 rk4 rk4_auto rk4_auto_midpoint);
@EXPORT_OK = qw(rk4_classical rk4_ralston arr2txt);
%EXPORT_TAGS = (ALL => [@EXPORT,@EXPORT_OK]);

sub new {
	my $arg1 = shift;
	my $class = ref($arg1) || $arg1; # can be used as class or instance method
	my $self  = {};   # ref to an empty hash
	bless $self, $class;
	$self->_initialise();
	return $self;
}

sub rk2 { my ($ynref, $dydtref, $t, $dt) = @_;
	if (ref $dydtref ne 'CODE') {
		warn "Math::RungeKutta::rk2: 2nd arg must be a subroutine ref\n";
		return ();
	}

	my $gamma = .75;  # Ralston's minimisation of error bounds
	my $alpha = 0.5/$gamma; my $beta = 1.0-$gamma;
	my $alphadt=$alpha*$dt; my $betadt=$beta*$dt; my $gammadt=$gamma*$dt;
	if (ref $ynref eq 'ARRAY') {
		my $ny = $#$ynref;
		my @ynp1; $#ynp1 = $ny;
		my @dydtn; $#dydtn = $ny;
		my @ynpalpha; $#ynpalpha = $ny;  # Gear calls this q
		my @dydtnpalpha; $#dydtnpalpha = $ny;
		@dydtn = &{$dydtref}($t, @$ynref);
		my $i; for ($i=$[; $i<=$ny; $i++) {
			$ynpalpha[$i] = ${$ynref}[$i] + $alphadt*$dydtn[$i];
		}
		@dydtnpalpha = &{$dydtref}($t+$alphadt, @ynpalpha);
		for ($i=$[; $i<=$ny; $i++) {
			$ynp1[$i]
			 = ${$ynref}[$i]+$betadt*$dydtn[$i]+$gammadt*$dydtnpalpha[$i];
		}
		return ($t+$dt, @ynp1);

	} elsif (ref $ynref eq 'HASH') {
		my %ynp1;
		my %ynpalpha;  # Gear calls this q
		my %dydtn = &{$dydtref}($t, %$ynref);
		foreach my $i (keys %$ynref) {
			$ynpalpha{$i} = ${$ynref}{$i} + $alphadt*$dydtn{$i};
		}
		my %dydtnpalpha = &{$dydtref}($t+$alphadt, %ynpalpha);
		foreach my $i (keys %$ynref) {
			$ynp1{$i}
			 = ${$ynref}{$i}+$betadt*$dydtn{$i}+$gammadt*$dydtnpalpha{$i};
		}
		return ($t+$dt, %ynp1);
	} else {
		warn "Math::RungeKutta::rk2: 1st arg must be an arrayref or hashref\n";
		return ();
	}
}
my @saved_k0; my %saved_k0; my $use_saved_k0 = 0;
sub rk4 { my ($ynref, $dydtref, $t, $dt) = @_;
	# The Runge-Kutta-Merson 5-function-evaluation 4th-order method
	# in the sine-cosine example, this seems to work as a 7th-order method !
	if (ref $dydtref ne 'CODE') {
		warn "Math::RungeKutta::rk4: 2nd arg must be a subroutine ref\n";
		return ();
	}
	if (ref $ynref eq 'ARRAY') {
		my $ny = $#$ynref; my $i;
		# @eta0 = @#ynref;
		my @k0; $#k0=$ny;
		if ($use_saved_k0) { @k0 = @saved_k0;
		} else { @k0 = &{$dydtref}($t, @$ynref);
		}
		for ($i=$[; $i<=$ny; $i++) { $k0[$i] *= $dt; }
		my @eta1; $#eta1=$ny;
		for ($i=$[; $i<=$ny; $i++) { $eta1[$i] = ${$ynref}[$i] + $k0[$i]/3.0; }
		my @k1; $#k1=$ny;
		@k1 = &{$dydtref}($t + $dt/3.0, @eta1);
		for ($i=$[; $i<=$ny; $i++) { $k1[$i] *= $dt; }
		my @eta2; $#eta2=$ny;
		my @k2; $#k2=$ny;
		for ($i=$[; $i<=$ny; $i++) {
			$eta2[$i] = ${$ynref}[$i] + ($k0[$i]+$k1[$i])/6.0;
		}
		@k2 = &{$dydtref}($t + $dt/3.0, @eta2);
		for ($i=$[; $i<=$ny; $i++) { $k2[$i] *= $dt; }
		my @eta3; $#eta3=$ny;
		for ($i=$[; $i<=$ny; $i++) {
			$eta3[$i] = ${$ynref}[$i] + ($k0[$i]+3.0*$k2[$i])*0.125;
		}
		my @k3; $#k3=$ny;
		@k3 = &{$dydtref}($t+0.5*$dt, @eta3);
		for ($i=$[; $i<=$ny; $i++) { $k3[$i] *= $dt; }
		my @eta4; $#eta4=$ny;
		for ($i=$[; $i<=$ny; $i++) {
			$eta4[$i] = ${$ynref}[$i] + ($k0[$i]-3.0*$k2[$i]+4.0*$k3[$i])*0.5;
		}
		my @k4; $#k4=$ny;
		@k4 = &{$dydtref}($t+$dt, @eta4);
		for ($i=$[; $i<=$ny; $i++) { $k4[$i] *= $dt; }
		my @ynp1; $#ynp1 = $ny;
		for ($i=$[; $i<=$ny; $i++) {
			$ynp1[$i] = ${$ynref}[$i] + ($k0[$i]+4.0*$k3[$i]+$k4[$i])/6.0;
		}
		# Merson's method for error estimation, see Gear p85, only works
		# if F is linear, ie F = Ay + bt, so that eta4 has no 4th-order
		# errors.  So in general step-doubling is the only way to do it.
		# Estimate error terms ...
		# if ($epsilon) {
		# 	my $errmax = 0; my $diff;
		# 	for ($i=$[; $i<=$ny; $i++) {
		# 		$diff = 0.2 * abs ($ynp1[$i] - $eta4[$i]);
		# 		if ($errmax < $diff) { $errmax = $diff; }
		# 	}
		# 	# print "errmax = $errmax\n"; # not much related to the actual error
		# }
		return ($t+$dt, @ynp1);

	} elsif (ref $ynref eq 'HASH') {
		my %k0;
		if ($use_saved_k0) { %k0 = %saved_k0;
		} else { %k0 = &{$dydtref}($t, %$ynref);
		}
		foreach my $i (keys(%$ynref)) { $k0{$i} *= $dt; }
		my %eta1;
		foreach my $i (keys(%$ynref)) { $eta1{$i} = ${$ynref}{$i}+$k0{$i}/3.0; }
		my %k1 = &{$dydtref}($t + $dt/3.0, %eta1);
		foreach my $i (keys(%$ynref)) { $k1{$i} *= $dt; }
		my %eta2;
		foreach my $i (keys(%$ynref)) {
			$eta2{$i} = ${$ynref}{$i} + ($k0{$i}+$k1{$i})/6.0;
		}
		my %k2 = &{$dydtref}($t + $dt/3.0, %eta2);
		foreach my $i (keys(%$ynref)) { $k2{$i} *= $dt; }
		my %eta3;
		foreach my $i (keys(%$ynref)) {
			$eta3{$i} = ${$ynref}{$i} + ($k0{$i}+3.0*$k2{$i})*0.125;
		}
		my %k3 = &{$dydtref}($t+0.5*$dt, %eta3);
		foreach my $i (keys(%$ynref)) { $k3{$i} *= $dt; }
		my %eta4;
		foreach my $i (keys(%$ynref)) {
			$eta4{$i} = ${$ynref}{$i} + ($k0{$i}-3.0*$k2{$i}+4.0*$k3{$i})*0.5;
		}
		my %k4 = &{$dydtref}($t+$dt, %eta4);
		foreach my $i (keys(%$ynref)) { $k4{$i} *= $dt; }
		my %ynp1;
		foreach my $i (keys(%$ynref)) {
			$ynp1{$i} = ${$ynref}{$i} + ($k0{$i}+4.0*$k3{$i}+$k4{$i})/6.0;
		}
		return ($t+$dt, %ynp1);
	} else {
		warn "Math::RungeKutta::rk4: 1st arg must be an arrayref or hashref\n";
		return ();
	}
}
my $t; my $halfdt;
my @y2; my %y2;  # need to be remembered for the midpoint
sub rk4_auto { my $ynref=shift; my $dydtref=shift; $t=shift;
	my ($dt, $arg4) = @_;
	if (ref $dydtref ne 'CODE') {
		warn "Math::RungeKutta::rk4_auto: 2nd arg must be a subroutine ref\n";
		return ();
	}
	if ($dt == 0) { $dt = 0.1; }
	if (ref $ynref eq 'ARRAY') {
		my @errors; my $epsilon;
		if (ref $arg4 eq 'ARRAY') {
			@errors = @$arg4; undef $epsilon;
		} else {
			$epsilon = abs $arg4; undef @errors;
			if (! $epsilon) { $epsilon = .0000001; }
		}
		my $ny = $#$ynref; my $i;
		my @y1; $#y1 = ny;
		$#y2 = ny; %y2 = ();
		my @y3; $#y3 = ny;
		$#saved_k0 = ny; @saved_k0 = &{$dydtref}($t, @$ynref);
		my $resizings = 0;
		my $highest_low_error = 0.1E-99; my $highest_low_dt = 0.0;
		my $lowest_high_error = 9.9E99;  my $lowest_high_dt = 9.9E99;
		while (1) {
			$halfdt = 0.5 * $dt; my $dummy;
			$use_saved_k0 = 1;
			($dummy, @y1) = &rk4($ynref, $dydtref, $t, $dt);
			($dummy, @y2) = &rk4($ynref, $dydtref, $t, $halfdt);
			$use_saved_k0 = 0;
			($dummy, @y3) = &rk4(\@y2, $dydtref, $t+$halfdt, $halfdt);
			my $relative_error;
			if ($epsilon) {
	 			my $errmax = 0; my $diff; my $ymax = 0;
	 			for ($i=$[; $i<=$ny; $i++) {
	 				$diff = abs ($y1[$i]-$y3[$i]);
	 				if ($errmax < $diff) { $errmax = $diff; }
	 				if ($ymax < abs ${$ynref}[$i]) {$ymax = abs ${$ynref}[$i];}
	 			}
				$relative_error = $errmax/($epsilon*$ymax);
			} elsif (@errors) {
				$relative_error = 0.0; my $diff;
	 			for ($i=$[; $i<=$ny; $i++) {
	 				$diff = abs ($y1[$i]-$y3[$i]) / abs $errors[$i];
	 				if ($relative_error < $diff) { $relative_error = $diff; }
	 			}
			} else { die
			 "RungeKutta::rk4_auto: \$epsilon & \@errors both undefined\n";
			}
			# Gear's correction assumes error is always in 5th-order terms :-(
			# $y1[$i] = (16.0*$y3{$i] - $y1[$i]) / 15.0;
			if ($relative_error < 0.60) {
				if ($dt > $highest_low_dt) {
					$highest_low_error = $relative_error;
					$highest_low_dt = $dt;
				}
			} elsif ($relative_error > 1.67) {
				if ($dt < $lowest_high_dt) {
					$lowest_high_error = $relative_error;
					$lowest_high_dt = $dt;
				}
			} else {
				last;
			}
			if ($lowest_high_dt<9.8E99 && $highest_low_dt>1.0E-99) { # interp
				my $denom = log ($lowest_high_error/$highest_low_error);
				if ($highest_low_dt==0.0||$highest_low_error==0.0||$denom==0.0){
					$dt = 0.5 * ($highest_low_dt+$lowest_high_dt);
				} else {
					$dt = $highest_low_dt * ( ($lowest_high_dt/$highest_low_dt)
				 	** ((log (1.0/$highest_low_error)) / $denom) );
				}
			} else {
				my $adjust = $relative_error**(-0.2); # hope error is 5th-order ...
				if (abs $adjust > 2.0) {
					$dt *= 2.0;  # prevent infinity if 4th-order is exact ...
				} else {
					$dt *= $adjust;
				}
			}
			$resizings++;
			if ($resizings>4 && $highest_low_dt>1.0E-99) {
				# hope a small step forward gets us out of this mess ...
				$dt = $highest_low_dt;  $halfdt = 0.5 * $dt;
				$use_saved_k0 = 1;
				($dummy, @y2) = &rk4($ynref, $dydtref, $t, $halfdt);
				$use_saved_k0 = 0;
				($dummy, @y3) = &rk4(\@y2, $dydtref, $t+$halfdt, $halfdt);
				last;
			}
		}
		return ($t+$dt, $dt, @y3);

	} elsif (ref $ynref eq 'HASH') {
		my %errors; my $epsilon;
		if (ref $arg4 eq 'HASH') {
			%errors = %$arg4; undef $epsilon;
		} else {
			$epsilon = abs $arg4; undef %errors;
			if (! $epsilon) { $epsilon = .0000001; }
		}
		my $i; my %y1; @y2 = (); my %y3;
		%saved_k0 = &{$dydtref}($t, %$ynref);
		my $resizings = 0;
		my $highest_low_error = 0.1E-99; my $highest_low_dt = 0.0;
		my $lowest_high_error = 9.9E99;  my $lowest_high_dt = 9.9E99;
		while (1) {
			$halfdt = 0.5 * $dt; my $dummy;
			$use_saved_k0 = 1;
			($dummy, %y1) = &rk4($ynref, $dydtref, $t, $dt);
			($dummy, %y2) = &rk4($ynref, $dydtref, $t, $halfdt);
			$use_saved_k0 = 0;
			($dummy, %y3) = &rk4(\%y2, $dydtref, $t+$halfdt, $halfdt);
			my $relative_error;
			if ($epsilon) {
	 			my $errmax = 0; my $diff; my $ymax = 0;
	 			foreach $i (keys(%$ynref)) {
	 				$diff = abs ($y1{$i}-$y3{$i});
	 				if ($errmax < $diff) { $errmax = $diff; }
	 				if ($ymax < abs ${$ynref}{$i}) {$ymax = abs ${$ynref}{$i};}
	 			}
				$relative_error = $errmax/($epsilon*$ymax);
			} elsif (%errors) {
				$relative_error = 0.0; my $diff;
	 			foreach $i (keys(%$ynref)) {
	 				$diff = abs ($y1{$i}-$y3{$i}) / abs $errors{$i};
	 				if ($relative_error < $diff) { $relative_error = $diff; }
	 			}
			} else { die
			 "RungeKutta::rk4_auto: \$epsilon & \%errors both undefined\n";
			}
			# Gear's correction assumes error is always in 5th-order terms :-(
			# $y1[$i] = (16.0*$y3{$i] - $y1[$i]) / 15.0;
			if ($relative_error < 0.60) {
				if ($dt > $highest_low_dt) {
					$highest_low_error = $relative_error;
					$highest_low_dt = $dt;
				}
			} elsif ($relative_error > 1.67) {
				if ($dt < $lowest_high_dt) {
					$lowest_high_error = $relative_error;
					$lowest_high_dt = $dt;
				}
			} else {
				last;
			}
			if ($lowest_high_dt<9.8E99 && $highest_low_dt>1.0E-99) { # interp
				my $denom = log ($lowest_high_error/$highest_low_error);
				if ($highest_low_dt==0.0||$highest_low_error==0.0||$denom==0.0){
					$dt = 0.5 * ($highest_low_dt+$lowest_high_dt);
				} else {
					$dt = $highest_low_dt * ( ($lowest_high_dt/$highest_low_dt)
				 	** ((log (1.0/$highest_low_error)) / $denom) );
				}
			} else {
				my $adjust = $relative_error**(-0.2); # hope error is 5th-order ...
				if (abs $adjust > 2.0) {
					$dt *= 2.0;  # prevent infinity if 4th-order is exact ...
				} else {
					$dt *= $adjust;
				}
			}
			$resizings++;
			if ($resizings>4 && $highest_low_dt>1.0E-99) {
				# hope a small step forward gets us out of this mess ...
				$dt = $highest_low_dt;  $halfdt = 0.5 * $dt;
				$use_saved_k0 = 1;
				($dummy, %y2) = &rk4($ynref, $dydtref, $t, $halfdt);
				$use_saved_k0 = 0;
				($dummy, %y3) = &rk4(\%y2, $dydtref, $t+$halfdt, $halfdt);
				last;
			}
		}
		return ($t+$dt, $dt, %y3);

	} else { die
		"Math::RungeKutta::rk4_auto: 1st arg must be arrayref or hashref\n";
		# return ();
	}
}
sub rk4_auto_midpoint {
	if (@y2) { return ($t+$halfdt, @y2); } else { return ($t+$halfdt, %y2); }
}

# ---------------------- EXPORT_OK routines ----------------------

sub rk4_ralston { my ($ynref, $dydtref, $t, $dt) = @_;
	if (ref $dydtref ne 'CODE') {
		warn "RungeKutta::rk4_ralston: 2nd arg must be a subroutine ref\n";
		return ();
	}
	# Ralston's minimisation of error bounds, see Gear p36
	if (ref $ynref eq 'ARRAY') {
		my $ny = $#$ynref; my $i;
		my $alpha1=0.4; my $alpha2 = 0.4557372542; # = .875 - .1875*(sqrt 5);
		my @k0; $#k0=$ny;
		@k0 = &{$dydtref}($t, @$ynref);
		for ($i=$[; $i<=$ny; $i++) { $k0[$i] *= $dt; }
		my @k1; $#k1=$ny;
		for ($i=$[; $i<=$ny; $i++) { $k1[$i] = ${$ynref}[$i] + 0.4*$k0[$i]; }
		@k1 = &{$dydtref}($t + $alpha1*$dt, @k1);
		for ($i=$[; $i<=$ny; $i++) { $k1[$i] *= $dt; }
		my @k2; $#k2=$ny;
		for ($i=$[; $i<=$ny; $i++) {
			$k2[$i] = ${$ynref}[$i] + 0.2969776*$k0[$i] + 0.15875966*$k1[$i];
		}
		@k2 = &{$dydtref}($t + $alpha2*$dt, @k2);
		for ($i=$[; $i<=$ny; $i++) { $k2[$i] *= $dt; }
		my @k3; $#k3=$ny;
		for ($i=$[; $i<=$ny; $i++) {
			$k3[$i] = ${$ynref}[$i] + 0.21810038*$k0[$i] - 3.0509647*$k1[$i]
		 	+ 3.83286432*$k2[$i];
		}
		@k3 = &{$dydtref}($t+$dt, @k3);
		for ($i=$[; $i<=$ny; $i++) { $k3[$i] *= $dt; }
		my @ynp1; $#ynp1 = $ny;
		for ($i=$[; $i<=$ny; $i++) {
			$ynp1[$i] = ${$ynref}[$i] + 0.17476028*$k0[$i]
		 	- 0.55148053*$k1[$i] + 1.20553547*$k2[$i] + 0.17118478*$k3[$i];
		}
		return ($t+$dt, @ynp1);

	} elsif (ref $ynref eq 'HASH') {
		my $i;
		my $alpha1=0.4; my $alpha2 = 0.4557372542; # = .875 - .1875*(sqrt 5);
		my %k0 = &{$dydtref}($t, %$ynref);
		foreach $i (keys(%$ynref)) { $k0{$i} *= $dt; }
		my %k1;
		foreach $i (keys(%$ynref)) { $k1{$i} = ${$ynref}{$i} + 0.4*$k0{$i}; }
		%k1 = &{$dydtref}($t + $alpha1*$dt, %k1);
		foreach $i (keys(%$ynref)) { $k1{$i} *= $dt; }
		my %k2;
		foreach $i (keys(%$ynref)) {
			$k2{$i} = ${$ynref}{$i} + 0.2969776*$k0{$i} + 0.15875966*$k1{$i};
		}
		%k2 = &{$dydtref}($t + $alpha2*$dt, %k2);
		foreach $i (keys(%$ynref)) { $k2{$i} *= $dt; }
		my %k3;
		foreach $i (keys(%$ynref)) {
			$k3{$i} = ${$ynref}{$i} + 0.21810038*$k0{$i} - 3.0509647*$k1{$i}
		 	+ 3.83286432*$k2{$i};
		}
		%k3 = &{$dydtref}($t+$dt, %k3);
		foreach $i (keys(%$ynref)) { $k3{$i} *= $dt; }
		my %ynp1;
		foreach $i (keys(%$ynref)) {
			$ynp1{$i} = ${$ynref}{$i} + 0.17476028*$k0{$i}
		 	- 0.55148053*$k1{$i} + 1.20553547*$k2{$i} + 0.17118478*$k3{$i};
		}
		return ($t+$dt, %ynp1);

	} else {
		warn "Math::RungeKutta::rk4_ralston: 1st arg must be arrayref or hashref\n";
		return ();
	}
}
sub rk4_classical { my ($ynref, $dydtref, $t, $dt) = @_;
	if (ref $dydtref ne 'CODE') {
		warn "RungeKutta::rk4_classical: 2nd arg must be subroutine ref\n";
		return ();
	}
	# The Classical 4th-order Runge-Kutta Method, see Gear p35
	if (ref $ynref eq 'ARRAY') {
		my $ny = $#$ynref; my $i;
		my @k0; $#k0=$ny;
		@k0 = &{$dydtref}($t, @$ynref);
		for ($i=$[; $i<=$ny; $i++) { $k0[$i] *= $dt; }
		my @eta1; $#eta1=$ny;
		for ($i=$[; $i<=$ny; $i++) { $eta1[$i] = ${$ynref}[$i] + 0.5*$k0[$i]; }
		my @k1; $#k1=$ny;
		@k1 = &{$dydtref}($t+0.5*$dt, @eta1);
		for ($i=$[; $i<=$ny; $i++) { $k1[$i] *= $dt; }
		my @eta2; $#eta2=$ny;
		for ($i=$[; $i<=$ny; $i++) { $eta2[$i] = ${$ynref}[$i] + 0.5*$k1[$i]; }
		my @k2; $#k2=$ny;
		@k2 = &{$dydtref}($t+0.5*$dt, @eta2);
		for ($i=$[; $i<=$ny; $i++) { $k2[$i] *= $dt; }
		my @eta3; $#eta3=$ny;
		for ($i=$[; $i<=$ny; $i++) { $eta3[$i] = ${$ynref}[$i] + $k2[$i]; }
		my @k3; $#k3=$ny;
		@k3 = &{$dydtref}($t+$dt, @eta3);
		for ($i=$[; $i<=$ny; $i++) { $k3[$i] *= $dt; }
		my @ynp1; $#ynp1 = $ny;
		for ($i=$[; $i<=$ny; $i++) {
			$ynp1[$i] = ${$ynref}[$i] +
		 	($k0[$i] + 2.0*$k1[$i] + 2.0*$k2[$i] + $k3[$i]) / 6.0;
		}
		return ($t+$dt, @ynp1);

	} elsif (ref $ynref eq 'HASH') {
		my %k0 = &{$dydtref}($t, %$ynref);
		foreach my $i (keys(%$ynref)) { $k0{$i} *= $dt; }
		my %eta1;
		foreach $i (keys(%$ynref)) { $eta1{$i} = ${$ynref}{$i} + 0.5*$k0{$i}; }
		my %k1 = &{$dydtref}($t+0.5*$dt, %eta1);
		foreach $i (keys(%$ynref)) { $k1{$i} *= $dt; }
		my %eta2;
		foreach $i (keys(%$ynref)) { $eta2{$i} = ${$ynref}{$i} + 0.5*$k1{$i}; }
		my %k2 = &{$dydtref}($t+0.5*$dt, %eta2);
		foreach $i (keys(%$ynref)) { $k2{$i} *= $dt; }
		my %eta3;
		foreach $i (keys(%$ynref)) { $eta3{$i} = ${$ynref}{$i} + $k2{$i}; }
		my %k3 = &{$dydtref}($t+$dt, %eta3);
		foreach $i (keys(%$ynref)) { $k3{$i} *= $dt; }
		my %ynp1;
		foreach $i (keys(%$ynref)) {
			$ynp1{$i} = ${$ynref}{$i} +
		 	($k0{$i} + 2.0*$k1{$i} + 2.0*$k2{$i} + $k3{$i}) / 6.0;
		}
		return ($t+$dt, %ynp1);

	} else {
		warn "Math::RungeKutta::rk4_classical: 1st arg must be arrayref or hashref\n";
		return ();
	}
}

# --------------------- infrastructure ----------------------

sub arr2txt { # neat printing of arrays for debug use
	my @txt; foreach (@_) { push @txt, sprintf('%g',$_); }
	return join (' ',@txt)."\n";
}
my $flag;
sub gaussn {   my $standdev = $_[$[];
	# returns normal distribution around 0.0 by the Box-Muller rules
	if (! $flag) {
		$a = sqrt(-2.0 * log(rand)); $b = 6.28318531 * rand;
		$flag = 1; return ($standdev * $a * sin($b));
	} else {
		$flag = 0; return ($standdev * $a * cos($b));
	}
}
1;

__END__