2.5.3 エネルギー保存則の確認

運動エネルギー、ポテンシャルエネルギーおよび全エネルギーを各時刻で 計算して、ファイルに出力する。

(1)

エネルギー保存則を確認するので、摩擦係数はゼロにしておく。 摩擦係数は入力ファイル中のgammaである。

(2)

出力ファイルを新しくopenする。ファイル名も指定する。 open文のところへもう一つのopen文を加える。

(3)

各時刻で運動エネルギー、ポテンシャルエネルギーおよび 全エネルギーを計算し、新しくopenしたファイルに出力する。 時間のdoループの中に

      ek(i)=v(i)*v(i)*am*0.5d0
      ev(i)=am*g*z(i)
      et(i)=ek(i)+ev(i)

としておき、これを新しく定義したファイルに出力すればよい。 配列ek,ev,etについては、配列宣言をしておかなければ 使えるようにならない。 座標の出力と同様にすればよいが、文番号(1から9999まで)については別の 整数を用いる。

(4)

出力したデータをグラフに描画する。 運動エネルギー、ポテンシャルエネルギー、全エネルギーの3つを同じ グラフに描画し、その変化について、考察する。摩擦が無いときには、 全エネルギーは保存されているだろうか。

(5)

エネルギー保存則の精度について全エネルギーの平均および 分散を計算し、画面上に出力する。

      if(abs(gamma).lt.1d-9)then
        etave=0d0
        etvar=0d0
        do i=1,nstep
          etave=etave+et(i)
        enddo
        etave=etave/nstep
        do i=1,nstep
          etvar=etvar+(et(i)-etave)**2
        enddo
        etvar=etvar/nstep
        write(06,950) etave,etvar
  950   format(' energy average :',e15.6,' variance :',e15.6)
      endif