振動子の数値計算プログラム例

下記のプログラムは、perlという簡易言語で書いたものです。 中身を解読しつつ、 自分の知っている言語(FortranやCやJava)に翻訳して動かしてみてください。


#!/usr/bin/perl -w
# ----------------------------------------------------------------
# 振動子の方程式を解く  (2000-11-16, 大信田 丈志)
#   dx/dt = v,   dv/dt = (1/m)F_spring(x) = -x,   m=1
# ----------------------------------------------------------------

# 刻み幅の設定
$dt = 0.05;

# 初期条件を読み込む
if( @ARGV == 2 ) {
    ($x0,$v0) = @ARGV;  $tmax = 6.5;
} elsif( @ARGV == 3 ) {
    ($x0,$v0,$tmax) = @ARGV;
} else {
    die "usage: $0 x0 v0 [tmax]\n";
}

# 数値計算
($x,$v) = ($x0,$v0);  $t = 0.0;

while( $t < $tmax ) {
    printf("%6.3f %f %f\n", $t, $x, $v);
    $xh   = $x + $v             *($dt/2);
    $vh   = $v + &F_spring($x)  *($dt/2);
    $xnew = $x + $vh            * $dt;
    $vnew = $v + &F_spring($xh) * $dt;
    ($x,$v) = ($xnew,$vnew);
    $t += $dt;
}

exit 0;

# ----------------------------------------------------------------
# バネの力
sub F_spring {
    my($x) = @_;
    return (-$x);   # [バネの力] = -kx,  k=1
}

# ----------------------------------------------------------------

[「常微分方程式」に戻る] [「常微分方程式のプログラム例」に戻る]