まちがいプログラム!

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


[力学演習のページに戻る] [数値計算のページに戻る]