22 lines
613 B
Perl
22 lines
613 B
Perl
sub runge_kutta {
|
|
my ($yp, $dt) = @_;
|
|
sub {
|
|
my ($t, $y) = @_;
|
|
my @dy = $dt * $yp->( $t , $y );
|
|
push @dy, $dt * $yp->( $t + $dt/2, $y + $dy[0]/2 );
|
|
push @dy, $dt * $yp->( $t + $dt/2, $y + $dy[1]/2 );
|
|
push @dy, $dt * $yp->( $t + $dt , $y + $dy[2] );
|
|
return $t + $dt, $y + ($dy[0] + 2*$dy[1] + 2*$dy[2] + $dy[3]) / 6;
|
|
}
|
|
}
|
|
|
|
my $RK = runge_kutta sub { $_[0] * sqrt $_[1] }, .1;
|
|
|
|
for(
|
|
my ($t, $y) = (0, 1);
|
|
sprintf("%.0f", $t) <= 10;
|
|
($t, $y) = $RK->($t, $y)
|
|
) {
|
|
printf "y(%2.0f) = %12f ± %e\n", $t, $y, abs($y - ($t**2 + 4)**2 / 16)
|
|
if sprintf("%.4f", $t) =~ /0000$/;
|
|
}
|