作成: 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
として求め直す]必要がある。