下記のプログラムは、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 } # ----------------------------------------------------------------