エネルギーの保存

　ルンゲ－クッタ法ではかなりの長時間，軌道を追跡することができたが，エネルギーが減少し続けてしまう． これは，いくら時間刻みを小さくしても改善しない．小さくしたなりに減少する．これをなんとかしたい．

Leap-frog法

蛙飛び法（Leap-frog scheme）を使ってみよう． 以下の微分方程式を計算するのに，

\[\begin{split}\dot x &= u, \dot u &= \Phi (x)\end{split}\]

次のようなアルゴリズムを用いる．

\[\begin{split}&u(t+\Delta t/2) = u(t-\Delta t/2) + \Phi (x) \Delta t \\ &x(t+\Delta t) = x(t) + u(t+\Delta t/2) \Delta t\end{split}\]

まず半歩先の速度を計算し，この速度を使って位置を更新する． tでの速度をが知りたいときは，以下の計算をする．

\[\begin{split}&u(t) = \{u(t-\Delta t/2) + u(t+\Delta t/2)\}/2\end{split}\]

また，最初のステップの計算を行う時には， \(u(-\Delta t/2)\) での値が必要である． これは以下のように計算する．

\[u(-\Delta t/2) = u(0) - \Phi (x) \Delta t/2\]

蛙飛び法はシンプレクィック解法と呼ばれるものの一種である． シンプレクィック解法は，長時間にわたってエネルギーの保存がよい（増えたり減ったりしない）， 時間反転対称性（時間を逆にたどると初期状態に戻る？）等の特徴を有する． 保存系の力学を解くのに優れている．2次精度の解法である．

Note

〔 計算物理学特論　分子動力学シミュレーション・数値積分法，能勢修一 〕ネットで入手

Note

〔 ハミルトン力学系のためのシンプレクティック数値積分法，吉田春夫 〕ネットで入手

!+++++++++++++++++++++++++++++++++++++++++++++++
module globals

double precision                delta_t

end module globals
!+++++++++++++++++++++++++++++++++++++++++++++++
module subfunc

implicit none

contains
    function fx(u,v) result(xx)

        double precision,intent(in)::    u,v
        double precision        xx

        xx = u

    end function fx
    function fu(x,y) result(xx)

        double precision,intent(in)::    x,y
        double precision        r,xx

        r = sqrt(x**2+y**2)
        xx = -x / r**3

    end function fu
    function fy(u,v) result(xx)

        double precision,intent(in)::    u,v
        double precision        xx

        xx = v

    end function fy
    function fv(x,y) result(xx)

        double precision,intent(in)::    x,y
        double precision        r,xx

        r = sqrt(x**2+y**2)
        xx = -y / r**3

    end function fv
end module subfunc
!+++++++++++++++++++++++++++++++++++++++++++++++
program celestial

use globals
use subfunc

implicit none

integer                     i,dout
integer,parameter::         n=10**8
integer,parameter::         noutmax=10000
double precision            tmax
double precision            t,x,y,u,v
double precision            uold,vold
double precision            u_a,u_b,v_a,v_b
double precision            E,E0,dE
double precision            start_time,finish_time

!-----CPU time
call CPU_TIME(start_time)

x = 1.0d0
y = 0.0d0
u = 0.0d0
v = 1.0d0
t = 0.0d0
tmax = 10.0d5

delta_t = (tmax-t)/dble(n)

open(unit=10,file='outputLF.dat',status='replace')
write(10,*)'t_lf, x_lf, y_lf, E_lf, dE_lf'

if(n .gt. noutmax) then
    dout = int(n/noutmax)
else
    dout = 1
end if

do i = 1, n-1

    if(i .eq. 1) then
        E0 = 0.5d0*(u**2+v**2) - 1.0d0/sqrt(x**2+y**2)
        u_b = u - 0.5d0*fu(x,y)*delta_t
        v_b = v - 0.5d0*fv(x,y)*delta_t
        uold = u
        vold = v
    end if

    u = uold
    v = vold

    E = 0.5d0*(u**2+v**2) - 1.0d0/sqrt(x**2+y**2)
    dE = abs((E - E0)/E0)

    if(mod(i-1,dout) .eq. 0) then
        write(10,*)t,',',x,',',y,',',E,',',dE
    end if

    u_a = u_b + fu(x,y)*delta_t
    v_a = v_b + fv(x,y)*delta_t
    x = x + u_a*delta_t
    y = y + v_a*delta_t
    uold = (u_a + u_b)/2.0d0
    vold = (v_a + v_b)/2.0d0
    u_b = u_a
    v_b = v_a

    t = t + delta_t

end do

close(10)

!-----Output CPU time
call CPU_TIME(finish_time)
write(*,*) 'CPU time:',finish_time - start_time

end program celestial
!+++++++++++++++++++++++++++++++++++++++++++++++

Note

時間が半歩ずれている位置と速度の計算のタイミングを合わせて，エネルギーの計算を行うなどを考えているうちに少々複雑になってしまった．もっとうまくかけるかも．

全エネルギーの時間変化でルンゲ－クッタ法と蛙飛び法の結果を比較しよう． 計算条件は共に，tmax = 1000000.0d0，n=100000000(時間刻み0.01)，さらに蛙飛びについては時間刻み0.001も行った．

時間の経過と共に，エネルギーが減少してしまうルンゲ－クッタ法の結果に対して， 蛙飛び法は，変動はあるものの初期値から離れない．この一定の変動幅は時間刻みを小さくすれば小さくなる． またn=100000000ステップ（時間刻み0.01）の場合で計算時間を比較するとルンゲ－クッタ法が23秒に対して，蛙飛び法は7秒程度である．

見方を変えて，エネルギーの初期値からの差を（初期値で規格化して）プロットした．縦軸は対数表示である． ルンゲ－クッタ法は初期値からどんどん離れていく一方で，蛙飛び法は初期値周辺の値を維持し続けるのがわかる． また，蛙飛び法について，時間刻みを小さくすることで差がグンと小さくなっていることも明確にわかる． 時間刻み0.001は0.01と比べて，ちょうど10 ^-3 程度小さい．0.01と0.1の間にも同様の差がある．

時間刻みの大きさと誤差の関係をみることで，数値積分が正しく行えているかどうかを確認することができる． 先に述べたように(>> 計算の精度 )例えば1次精度のオイラー法では \(O((\Delta t)^2)\) を計算しない．この打ち切り誤差は時間刻みを10分の1にすれば，100分の1になるはずである． 蛙飛び法は2次精度なので， \(O((\Delta t)^3)\) を計算しない．時間刻み10分の1で誤差は10 ^-3 になるはずである． 時間刻みと誤差の間にこの通りの関係があるかどうかが，プログラムが正しくかけているかどうかの判断の一つになる．

つぎのような量を計算してもよい．

\[\frac{1}{N}\sum_{i}^N{|E_{i+1}-E_i|}\]

1ステップを更新する毎のエネルギー変化の平均値である． 蛙飛び法での計算において，上記の量を計算したところ \(\Delta t\) が0.1で1.58*10 ^-4 ， 0.01で1.59*10 ^-7 ，0.001で1.59*10 ^-10 であった．時間刻み10分の1で誤差は10 ^-3 になる．

Note

〔 分子動力学法入門，能勢修一 〕ネットで入手

速度ベルレ法

速度ベルレ法（Velocity Verlet scheme）も試してみよう． 以下の微分方程式を計算するのに，

\[\begin{split}\dot x &= u, \dot u &= \Phi (x)\end{split}\]

次のようなアルゴリズムを用いる．

\[\begin{split}&x(t+\Delta t) = x(t) + u(t) \Delta t + \Phi (x) \frac{(\Delta t)^2}{2}\\ &u(t+\Delta t) = u(t) + \{ \Phi (x(t)) + \Phi (x(t+\Delta t)) \}\frac{\Delta t}{2}\end{split}\]

速度ベルレ法もシンプレクィック解法である．？○？次精度である． 蛙飛び法と違って，位置と速度の更新タイミングが揃っている． プログラムの主要部分のみ以下に示す．

fu_temp = fu(x,y)
fv_temp = fv(x,y)
x = x + fx(u,v)*delta_t + 0.5d0*fu_temp*delta_t**2
y = y + fy(u,v)*delta_t + 0.5d0*fv_temp*delta_t**2
u = u + 0.5d0*(fu_temp + fu(x,y))*delta_t
v = v + 0.5d0*(fv_temp + fv(x,y))*delta_t

t = t + delta_t

初期値からのエネルギーの差の図に，時間刻み0.01と0.001の速度ベルレ法の結果を追加した．

同じ時間刻みで比較すると，一番良い結果が得られている． 1ステップを更新する毎のエネルギー変化の平均値は， \(\Delta t\) が0.1で3.91*10 ^-7 ， 0.01で3.98*10 ^-12 ，0.001で1.05*10 ^-16 であった． 時間刻み10分の1に対する誤差の差の規則性が0.01から0.001で維持されない．

Note

図からも0.001の速度ベルレでその他とは違う傾向が分かる．誤差自体はとても小さいが，特徴の一つであったはずの初期値からのずれが時間変動してしまっている．なぜだろう？打切り誤差に加えて丸め誤差の影響が現れたか？0.0001は0.001とほぼ同様の傾向．誤差が全く小さくならない．「double precisionは倍精度実数型の宣言．倍精度とは約15桁の精度で表現される実数である」

Note

蛙飛び法や速度ベルレ法など（陽的なシンプレクティック積分法）は，ハミルトニアンが位置の関数と速度の関数に分解できる時に成立する．簡単に言えば，全エネルギーが運動エネルギー（速度の関数）と位置エネルギー（位置の関数）の和で記述できるような場合である．