作成: 2000-11-15
どうも間違った数値計算プログラムが流布しているようだ。 たとえば、以下に示すプログラムは、2000年10月の課題にある 惑星の運動の問題を解こうとしたものらしいのだが…?
program rikigaku2
parameter(n=300000)
integer i
real h,fl,xah,vxh,xh,f2,yah,vyh,yh,r
real t(n),x(n),vx(n),y(n),vy(n)
h=0.0001
vx(1)=0.000
x(1)=0.500
vy(1)=1.630
y(1)=0.000
t(1)=0.000
do i=1,45000
r=sqrt(x(i)**2+y(i)**2)
f1=-x(i)/r**3
xh=x(i)+vx(i)*(0.3)*h
vxh=vx(i)+f1*(0.3)*h
xah=-xh/r**3
vx(i+1)=vx(i)+h*xah
x(i+1)=x(i)+h*vxh
f2=-y(i)/r**3
yh=y(i)+vy(i)*(0.3)*h
vyh=vy(i)+f2*(0.3)*h
yah=-yh/r**3
vy(i+1)=vy(i)+h*yah
y(i+1)=y(i)+h*vyh
t(i+1)=t(i)+h
enddo
do i=1,45000
write(*,*)t(i),x(i),vx(i),y(i),vy(i)
enddo
stop
end
このプログラムは、
そもそも数値計算の式が間違っている。
途中に0.3などとあるのはどういうことか?
うわさでは「大信田先生が0.3にしろと言った」ということに
なっているようだが、私はそんなことを言ったおぼえはない。
どうも話が間違って伝わっているようである。
まずは改良Euler法の説明をきちんと読んでほしい。
そして、この方法を 惑星の問題に正しくあてはめる必要がある。
具体的な方法については、間違いのないように、よく考えること!
正しくプログラムを組めば、
h=0.01程度でも ぴったりもとの位置に戻ってくるはずだ。
上記の例には、(「0.3」のほかにも)問題点がたくさんあるが、 そのなかからひとつだけ指摘しておく:
xah=-xh/r**3
これはおそらく、中間時刻における加速度を計算しているのであろう。
だが、rの中身の値が古いままなのはどういうわけか?
正しくは、
中間時刻における位置座標(xh, yh)を求めた時点で、
rも計算しなおす[またはrhとして求め直す]必要がある。