======================
計算の精度
======================

.. toctree::
   :maxdepth: 1


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

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

そもそも
=============

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

.. math::
    f(t+\Delta t)= f(t) + \Phi (t) \Delta t= f(t) + \frac{df}{dt} \Delta t
    :label: euler_scheme

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

.. index:: Tayer expansion

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

.. math::
    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 + ...
    :label: tayler expansion

| である． :eq:`euler_scheme` と :eq:`tayler expansion` を比較すると， :eq:`tayler expansion` の右辺第３項以降が :eq:`euler_scheme` にはないことがわかる．
| この右辺第３項以降を :math:`O((\Delta t)^2)` と書くとテイラー展開は以下のようになる．

.. math::
    f(t+\Delta t) = f(t) + \frac{df}{dt} \Delta t + O((\Delta t)^2)

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

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

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

差を縮める１
===================

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


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

.. figure:: ../img/308XvsYn100.PNG
   :scale: 100 %
   :align: center
   :alt: X vs Y

   Fig. :math:`\Delta t` dependence on trajectory

.. figure:: ../img/309TotalEn100.PNG
   :scale: 100 %
   :align: center
   :alt: Total energy

   Fig. :math:`\Delta t` dependence on total energy

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

.. index:: 
    single: scheme; Runge-Kutta

差を縮める２
===================

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

.. figure:: ../img/310XvsY_RK.PNG
   :scale: 100 %
   :align: center
   :alt: X vs Y

   Fig. Trajectory　by Runge-Kutta

.. figure:: ../img/311TotalE_RK.PNG
   :scale: 100 %
   :align: center
   :alt: Total energy

   Fig. :Total energy by Runge-Kutta

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

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

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

.. index:: file open

.. code-block:: fortran

    !+++++++++++++++++++++++++++++++++++++++++++++++
    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
    !+++++++++++++++++++++++++++++++++++++++++++++++