計算の精度

　得られたオイラー法による計算結果は軌道が解析解から少しずれ，全エネルギーが徐々に増加した．

なぜだろう．またこの解析解との差を縮めるためには何ができるだろうか．

そもそも

オイラー陽解法では以下の様な計算を行っている．

(1)\[f(t+\Delta t)= f(t) + \Phi (t) \Delta t= f(t) + \frac{df}{dt} \Delta t\]\[ [ \dot f (= \frac{df}{dt}) =\Phi (t) ]\]

ここでテイラー展開（Tayler expansion）を思い出すと，

(2)\[f(t+\Delta t) = f(t) + \frac{df}{dt} \Delta t + \frac{d^2 f}{dt^2} (\Delta t)^2 + \frac{d^3 f}{dt^3} (\Delta t)^3 + ... + \frac{d^n f}{dt^n} (\Delta t)^n + ...\]

である． (1) と (2) を比較すると， (2) の右辺第３項以降が (1) にはないことがわかる．

この右辺第３項以降を \(O((\Delta t)^2)\) と書くとテイラー展開は以下のようになる．

\[f(t+\Delta t) = f(t) + \frac{df}{dt} \Delta t + O((\Delta t)^2)\]

ということで，オイラー陽解法の計算手順は \(O((\Delta t)^2)\) の分だけテイラー展開と差があるということである．

Note

テイラー展開は，関数 \(f(t)\) が無限回微分可能で無限回計算を続けることができるならば，厳密に等号が成り立つと言ってよい（と思う）（例外はある）．オイラー陽解法はテイラー展開の計算を途中で（第２項目までで）打ち切っている．この打ち切りによる誤差を打切り誤差という．２次以上の項 \(O((\Delta t)^2)\) を無視しているので，１次精度の解法と呼ばれる．

Note

\(\Delta t\) を極限まで小さくできるならば，オイラー法陽解法もテイラー展開も計算結果は同じになる．しかし，コンピュータを用いた計算では \(\Delta t\) は値を持っていなくてはいけない．“極限まで小さい”値はもてない．

差を縮める１

\(O((\Delta t)^2)\) は \(\Delta t\) が小さければ小さいほど小さくなる． 先のプログラムコートでは，全計算時間をn=10分割して \(\Delta t\) を決めていたので，nをより大きくすれば \(\Delta t\) が小さくできる．

n=10を100に変えて計算した時の軌道と全エネルギーの変化は以下のようになる．

Fig. \(\Delta t\) dependence on trajectory

Fig. \(\Delta t\) dependence on total energy

どちらの図からも，n=100にすることで解析解に一気に近づいていることが分かる． ただし，nを大きくするということは，その分だけdoループによる反復計算の回数が増すということであり，計算に要する時間が増す． nを10倍（ \(\Delta t\) を10分の1）にすれば計算時間は10倍である．

差を縮める２

オイラー陽解法の精度に不満があるならば，計算方法（scheme）自体を変更するという選択肢もある． 良く知られている方法に4次のルンゲ-クッタ法（4th order Runge-Kutta scheme）がある． 以下は，この方法を用いてボール投げ問題を計算した結果である．オイラー法の結果と比較している． \(\Delta t\) はオイラー法のn=10の場合と同じである．

Fig. Trajectory　by Runge-Kutta

Fig. :Total energy by Runge-Kutta

この図の範囲では，解析解に完全に一致しているように見える．プログラムコードは若干複雑になり， 計算コストもそれなりに増すが，オイラー法でn=100とした場合よりもはるかによい一致である．

Note

ルンゲ－クッタ法の詳細は省略（いつかそのうち）．ルンゲ－クッタ法は４次精度である．

一応，コードも載せておく，．さきほどコードが若干複雑になると言った割には大幅な変更を受けいるように見える． モジュールや関数を新たに使っているためである．コード内には色々無駄がある．関数内の時間tは使っていない．

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

double precision,parameter::    g = 9.8d0
double precision,parameter::    pi = acos(-1.0d0)
double precision                delta_t

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

use globals
implicit none

contains
    function fx(t,x,u) result(xx)

        double precision,intent(in)::    t,x,u
        double precision        xx

        xx = u

    end function fx
    function fu(t,x,u) result(xx)

        double precision,intent(in)::    t,x,u
        double precision        xx

        xx = 0.0d0

    end function fu
    function fy(t,x,u) result(xx)

        double precision,intent(in)::    t,x,u
        double precision        xx

        xx = u

    end function fy
    function fv(t,x,u) result(xx)

        double precision,intent(in)::    t,x,u
        double precision        xx

        xx = -g

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

use globals
use subfunc

implicit none

integer                     i
integer,parameter::         n=10
double precision            v0,theta,tmax
double precision            t(n),x(n),y(n),u(n),v(n)
double precision            xa(n),ya(n),ua(n),va(n)
double precision            E(n),Ea(n)
double precision            kx1,kx2,kx3,kx4,ku1,ku2,ku3,ku4
double precision            ky1,ky2,ky3,ky4,kv1,kv2,kv3,kv4
double precision            start_time,finish_time

v0 = 10.0d0
theta = 30.0d0
theta = theta * pi / 180.0d0    !degree->rad
x(1) = 0.0d0
y(1) = 0.0d0
u(1) = v0 * cos(theta)
v(1) = v0 * sin(theta)
t(1) = 0.0d0
tmax = 1.0d0

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

do i = 1, n-1
    kx1 = fx(t(i),x(i),u(i))*delta_t
    ku1 = fu(t(i),x(i),u(i))*delta_t
    kx2 = fx(t(i)+0.5d0*delta_t,x(i)+0.5d0*kx1,u(i)+0.5d0*ku1)*delta_t
    ku2 = fu(t(i)+0.5d0*delta_t,x(i)+0.5d0*kx1,u(i)+0.5d0*ku1)*delta_t
    kx3 = fx(t(i)+0.5d0*delta_t,x(i)+0.5d0*kx2,u(i)+0.5d0*ku2)*delta_t
    ku3 = fu(t(i)+0.5d0*delta_t,x(i)+0.5d0*kx2,u(i)+0.5d0*ku2)*delta_t
    kx4 = fx(t(i)+delta_t,x(i)+kx3,u(i)+ku3)*delta_t
    ku4 = fu(t(i)+delta_t,x(i)+kx3,u(i)+ku3)*delta_t

    ky1 = fy(t(i),y(i),v(i))*delta_t
    kv1 = fv(t(i),y(i),v(i))*delta_t
    ky2 = fy(t(i)+0.5d0*delta_t,y(i)+0.5d0*ky1,v(i)+0.5d0*kv1)*delta_t
    kv2 = fv(t(i)+0.5d0*delta_t,y(i)+0.5d0*ky1,v(i)+0.5d0*kv1)*delta_t
    ky3 = fy(t(i)+0.5d0*delta_t,y(i)+0.5d0*ky2,v(i)+0.5d0*kv2)*delta_t
    kv3 = fv(t(i)+0.5d0*delta_t,y(i)+0.5d0*ky2,v(i)+0.5d0*kv2)*delta_t
    ky4 = fy(t(i)+delta_t,y(i)+ky3,v(i)+kv3)*delta_t
    kv4 = fv(t(i)+delta_t,y(i)+ky3,v(i)+kv3)*delta_t

    x(i+1) = x(i) + (kx1+2.0d0*kx2+2.0d0*kx3+kx4)/6.0d0
    u(i+1) = u(i) + (ku1+2.0d0*ku2+2.0d0*ku3+ku4)/6.0d0
    y(i+1) = y(i) + (ky1+2.0d0*ky2+2.0d0*ky3+ky4)/6.0d0
    v(i+1) = v(i) + (kv1+2.0d0*kv2+2.0d0*kv3+kv4)/6.0d0
    t(i+1) = t(i) + delta_t
end do

xa= x(1) + u(1)*t
ya = y(1) + v(1)*t - 0.5d0*g*t**2
ua = u(1)
va = v(1) - g*t
E = 0.5d0*(u**2 + v**2) + g*y
Ea = 0.5d0*(ua**2 + va**2) + g*ya

open(unit=10,file='outputRK.dat',status='replace')
write(10,*)'t, xrk, yrk, xa, ya, Erk, Ea'
do i = 1,n
    write(10,*)t(i),',',x(i),',',y(i),',',xa(i),',',ya(i),',',E(i),',',Ea(i)
end do
close(10)

end program ball_throw
!+++++++++++++++++++++++++++++++++++++++++++++++