2.6.1 関数副プログラムの導入

ルンゲ-クッタ法は、計算時間よりも計算精度を重要視するときに使用される。 計算するべき粒子に働く力を、1ステップ当たり4回計算しなければならないので プログラムも当然複雑になる。そこでプログラミングの見通しを良くして、 プログラミングの誤りを抑えるために、力の計算を関数副プログラム内にて行う。

元のプログラムでは、力は、

      f(i)=-am*g

と計算される。これは、時刻、位置座標、速度座標に依存していないが、 一般的には依存している。そこで、ルンゲ-クッタ法で必要となる加速度(力)と 速度の計算に関数を導入する。 (講義では、関数$f_{1}(t,x,v)$と$f_{2}(t,x,v)$として登場していたはずだ。) ここでは、計算式に質量やばね定数あるいはその他に関係 する定数が存在しているので、その全てを関数の引数にして関数を導入すること にしよう。 例えば、座標の更新を行うときに、 関数副プログラムfaとfvを導入して、主プログラムのでは、 methodでif文分岐するところに

      elseif(method.eq.1)then
      do i=1,nstep
        t=dt*(i-1)
        t1=t+dt*0.5d0
        vi=v(i)
        zi=z(i)
        xk1=dt*fa(t,vi,zi,am,g,ak0,hml,gamma)
        xl1=dt*fv(t,vi,zi,am,g,ak0,hml,gamma)
        xk2=dt*fa(t1,vi+xk1*0.5d0,zi+xl1*0.5d0,am,g,ak0,hml,gamma)
        xl2=   .
        xk3=   .
        xl3=   .
        xk4=   .
        xl4=   .
        v(i+1)=vi+(xk1+xk2*2d0+xk3*2d0+xk4)/6d0
        z(i+1)=zi+(xl1+xl2*2d0+xl3*2d0+xl4)/6d0
      enddo

とような感じでプログラミングしておけばよい。 また関数副プログラムfaとfv を主プログラム(fcode-bungy.f)の後に付けておく。 関数副プログラムfaは、例えば、

      double precision function fa(t,v,z,am,g,ak0,hml,gamma)
c
c     copyright Tatsuki Oda (2009.05.01)
c
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
c
      if(z.gt.hml)then
        ak=0d0
      else
        ak=ak0
      endif
      fa=-g-(ak*(z-hml)+gamma*v)/am
c
      return
      end

(関数)副プログラムにはreturnを入れることができて、 これに制御が移るとその(関数)副プログラムが終了して 値をfaに入れて返す。関数副プログラムfvについても 作成してプログラムの後に付けておく。

プログラムfcode-bungy.fを変更したので、コンパイルしておく。 入力データのmethodの指定値を1にして実行することにより ルンゲ-クッタ法で計算できるようになった。

method=0と同じ機能を持たせるにはev(i),ev(i),et(i) が計算されるように必要な部分を追加する必要がある。例えば、ev(i) の計算では、ak(ばね定数)が適切に設定されるように計算式を 追加する必要がある。

結果の確認を行う。methodの値を0と1とでそれぞれ1回ずつ 実行することにより、計算結果のアルゴリズムの違いによる比較を することができる。 その際に、1つ目を実行したときに、出力ファイルの名前をmvコマンド にて変更してけば、gnuplot上でデータの比較が可能となろう。