===========================
オイラー法でボール投げ
===========================

.. toctree::
   :maxdepth: 1

　ボールの運動を追跡する．水平方向からの角度 :math:`\theta` ，初速度 :math:`v_0` で投げたボールの軌跡が知りたい．

.. index:: equation of motion

運動方程式
=====================

| 運動について何か知りたければ，まずはとにかく運動方程式である．
| 水平方向座標を :math:`x` ，高さ方向座標を :math:`y` とし，時間 :math:`t` とする．ボールの質量が :math:`m` で重力加速度を :math:`g` のとき，
| 運動方程式は以下のように書ける．

.. math:: 
    m\frac{d^2 x}{dt^2} &= 0 \\
    m\frac{d^2 y}{dt^2} &= -mg
    :label: eq_motion

| もちろんこの方程式は解析的に解を求める（関数を使って解を記述する）ことができる．時間 :math:`t` で積分すればよい．

.. math::
    x &= x_0 + v_0\cos{\theta}\ t \\
    y &= y_0 + v_0\sin{\theta}\ t - \frac{1}{2}gt^2
    :label: analytical solution

| 時間 :math:`t` を指定すれば，既知の情報〔初期位置（:math:`x_0`,　:math:`y_0`）,初速度（:math:`v_0`, :math:`\theta`）, :math:`g`〕を使って位置を計算できる．速度は上式を時間で微分すればよい．


.. index:: 
    single: scheme; Euler's explicit

数値積分（オイラー法）
===============================

解析解が分かっているのだから本来必要ないが，数値的にも運動方程式を解いてみよう．

以下のような微分方程式について

.. math::
    \dot f = \Phi (t)
    :label: differential

左辺の微分演算は，

.. math:: 
    \dot f = \frac{df}{dt} = \lim\limits_{\epsilon \to 0} \frac{f(t+\epsilon)-f(t)}{\epsilon} 

であるが，これを

.. math::
    \lim\limits_{\epsilon \to 0} \frac{f(t+\epsilon)-f(t)}{\epsilon} \approx \frac{f(t+\Delta t)-f(t)}{\Delta t}

| とする． :math:`\Delta t` が小さいのであれば，このように近似してもよいだろう．
| すると，元の式は，

.. math::
    \dot f = \frac{f(t+\Delta t)-f(t)}{\Delta t} = \Phi (t)

となる．これを変形すると，

.. math::
    f(t+\Delta t) = f(t) + \Phi (t) \Delta t
    :label: euler_scheme

| この式は， :math:`f(t)`,  :math:`\Phi (t)`, :math:`\Delta t` が分かっていれば， :math:`f(t+\Delta t)` が計算できることを意味する．
| つまり，時刻 :math:`t` での情報 :math:`f(t)` , :math:`\Phi (t)` と時間刻み :math:`\Delta t` により，時刻 :math:`t + \Delta t` での情報  :math:`f(t + \Delta t)` が得られる．
| この計算方法をオイラー陽解法（Euler's explicit scheme）という．

.. note::
    :math:`t + \Delta t` での情報を得るのに， :math:`t` での情報のみを用いて計算する方法を陽解法という．そうではない，:math:`t + \Delta t` での情報も使って計算する方法を陰解法（implicit scheme）という．

では，計算してみよう
===============================

計算手順
-----------------

式 :eq:`differential` は一階の微分方程式で，式 :eq:`eq_motion` は二階である．
このままではオイラー陽解法を適用することができない．
そこで新たな変数 :math:`u = \dot x = \frac{dx}{dt}, v = \dot y` を導入しよう．

.. note::
    :math:`\dot x = \frac{dx}{dt}` のように，変数の上にドット「 :math:`\dot o` 」をつけるときは，時間微分を意味する，ことにする．

すると，式 :eq:`eq_motion` の二つの式は以下のように書ける．

.. math::
    \dot x &= u, \dot u &= 0 \\
    \dot y &= v, \dot v &= -g

式の数は増えてしまったが，いずれも一階の微分方程式でありオイラー法 :eq:`euler_scheme` が適用できる．

.. math::
    u(t+\Delta t) &= u(t) + (0) \Delta t = u(t)\\
    x(t+\Delta t) &= x(t) + u(t) \Delta t \\
    v(t+\Delta t) &= v(t) + (-g) \Delta t = v(t) -g \Delta t\\
    y(t+\Delta t) &= y(t) + v(t) \Delta t
    :label: euler_throw

時刻 :math:`t` での情報 :math:`u(t),x(t),v(t),y(t)` と時間刻み :math:`\Delta t` により，
時刻 :math:`t + \Delta t` での情報  :math:`u(t + \Delta t),x(t + \Delta t)`, :math:`v(t + \Delta t),y(t + \Delta t)` が得られる．

さて時刻 :math:`t` での情報 :math:`u(t),x(t),v(t),y(t)` とは？と思うかもしれない．しかし，問題設定より「角度 :math:`\theta` ，初速度 :math:`v_0` で投げる」ということを知っている． 
これより :math:`t=0` で， :math:`u(0)=v_0 \cos \theta` , :math:`v(0)=v_0 \sin \theta` が分かる．
また時刻 :math:`t=0` でのボールの位置を座標の原点とすれば， :math:`x(0)=0, y(0)=0` である．
つまり，  :math:`t=0` での情報 :math:`u(0),x(0)`, :math:`v(0),y(0)` が決まっているので，
これらの情報から :math:`t=0+\Delta t=\Delta t` の情報 :math:`u(\Delta t),x(\Delta t),v(\Delta t),y(\Delta t)` が計算できる． 
さらに :math:`t=\Delta t` での情報を使えば  :math:`t=\Delta t+\Delta t=2\Delta t` の情報が計算できる． 
この計算を順次n回繰り返せば，  :math:`t=0` での情報（ **初期条件** という）から出発して  :math:`t=n\Delta t` までの情報が計算できる．
これはつまり任意の時間での位置と速度の情報が得られるということに他ならない．（数値的に）運動方程式が解けたのである． 

以上の手順は，簡単な足し算・引き算と掛け算の演算のみで計算できる．もちろん手計算でやってもよい．．．

Fortranで計算しよう
-----------------------------

| コンピュータは単純な計算を繰り返し行うことが得意である．
| Fortranを使って計算手順を記述し，コンピュータで計算を実行してみよう．

.. code-block:: fortran
    :linenos:
    :emphasize-lines: 27-33
    
    !+++++++++++++++++++++++++++++++++++++++++++++++
    program ball_throw
    !+++++++++++++++++++++++++++++++++++++++++++++++

    implicit none

    integer                     i
    integer,parameter::         n=10
    double precision,parameter::    g = 9.8d0
    double precision,parameter::    pi = acos(-1.0d0)        
    double precision            v0,theta,tmax
    double precision            delta_t
    double precision            t(n),x(n),y(n),u(n),v(n)

    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
        u(i+1) = u(i)
        v(i+1) = v(i) - g*delta_t
        x(i+1) = x(i) + u(i)*delta_t
        y(i+1) = y(i) + v(i)*delta_t
        t(i+1) = t(i) + delta_t
    end do

    write(*,*)'t, x, y'
    do i = 1,n
        write(*,*)t(i),',',x(i),',',y(i)
    end do 

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

オイラー陽解法の主要部分はハイライトされている27-33行目である．しかし，他にも色々な手続きを書く必要がある．
pythonになれていると，面倒くさいと思ってしまうが仕方がない．

.. code-block:: fortran

    !***********

「!」に続く文字はプログラム内で無視される．コメントアウトである．プログラム内に説明書きを加えるたり，
区切りを入れたい時に使う．

.. code-block:: fortran

    program ball_throw


    end program ball_throw

「program」ではじまり，「end program」で終わる．プログラムの名前を付ける．

.. index:: implicit none

.. code-block:: fortran

    implicit none

暗黙の型宣言を行わない，の意味．暗黙の型宣言が有効な場合，i-nの頭文字で始まる変数を整数型，
それ以外が実数型になる．間違いの元になるので，変数はすべて型を宣言してから使う．

.. code-block:: fortran

    integer                     i
    integer,parameter::         n=10

プログラム内で使用する変数の型宣言．iは整数型（integer）です．としている．
「parameter」は，変数が定数であることを宣言する．プログラム内で間違えて値を変えようとすると警告してくれる．
「::」は変数の値を予め指定する場合に必要． つまり「integer,parameter::   n=10」は，変数nは値10の整数で定数であると宣言している．

.. index:: doble precision

.. code-block:: fortran

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

double precisionは倍精度実数型の宣言．倍精度とは約15桁の精度で表現される実数である．
単精度（精度6桁）も使えるが，どうしてもという時以外は倍精度でよい
（単精度を使えば変数を使うのに必要な記憶メモリの量は減らすことができる）．
gという変数に倍精度実数で9.8を意味する9.8d0を代入している．
9.8d0=9.8*10\ :sup:`0` \=9.8である．gは重量加速度であり，変化する値ではないのでparameterがついている．
次のpi = acos(-1.0d0)では，piという変数にacos(-1.0d0)で計算した値を入れている．このように計算した結果を使うこともできる．
acos(x)はcos\ :sup:`-1` \(x)を計算する関数である．
つまりpiにcos\ :sup:`-1` \(-1)の計算結果を入れている．円周率 :math:`\pi` を計算して求めている．
もちろんpiに値をいれるのにpi=3.141592....d0と数値を使ってもよい．

.. code-block:: fortran

    double precision            t(n),x(n),y(n),u(n),v(n)

倍精度配列の確保である．たとえばt(n)は，先にn=10としてあるので，t(1)~t(10)の10個の実数の確保を意味する．

.. code-block:: fortran

    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

計算に必要な初期条件を入力している．
v0は初速度を10m/s，thetaは角度30度，x(1)はx方向の初期位置，u(1)はx方向の初期速度である．
t(1)は計算開始時刻．tmaxは計算終了時刻を決めている．
cosやsin関数（先のacosも）は引数をラジアン単位にする必要があるので，17行目で度単位で入力したthetaをラジアン単位に変換している．

.. code-block:: fortran

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

tmax-t(1)は計算終了時刻と計算開始時刻の差なので，全計算時間である．これをdble(n)で割り算している．
dbleは引数を倍精度実数に変換する関数である．元々nは整数の10であるが，dble(10)=1.0d0，倍精度実数の10.0に変換される．
精度に差のある型の変数を用いた演算を行う場合は，強い型の精度が採用される．実数は整数よりも強いので，ここでは倍精度になる．
しかしここでは安全のために（というか癖だが），整数を明示的に倍精度実数に変換して割り算を行った．
ということで，delta_tは :math:`\Delta t` であり，(1.0-0.0)/10.0=0.1である．
nを大きくすればdelta_tを小さくできる．

.. index:: 
    pair: fortran; loop

.. code-block:: fortran

    do i = 1, n-1
        u(i+1) = u(i)
        v(i+1) = v(i) - g*delta_t
        x(i+1) = x(i) + u(i)*delta_t
        y(i+1) = y(i) + v(i)*delta_t
        t(i+1) = t(i) + delta_t
    end do

さて，やっといよいよオイラー陽解法の主要部分である．
「do i = 1, n-1」は反復計算の開始を宣言する．
i が1からn-1になるまで，「end do」までに書かれている処理を繰り返す．
ここでは，n=10でなので，i=1,2,...,8,9まで，9回の反復計算を行う．
反復計算の中身は．．．じっと見ていればわかると思うが， :eq:`euler_throw` の計算である．
「=」の右側（例えばv(i) - g*delta_t）に書かれている計算結果を「=」の左側（v(i+1)）の変数に代入している．
uからtまで一通りの計算が済んだら，1大きくなったiで再度計算を繰りかえす．
i=1のとき，v(1)（これは初期条件で与えらえている）を使ってv(1+1=2)を計算する．
i=2のとき，v(2)（i=1で計算してある）を使ってv(2+1=3)を計算する．

.. code-block:: fortran

    write(*,*)'t, x, y'
    do i = 1,n
        write(*,*)t(i),',',x(i),',',y(i)
    end do 

最後に，計算結果を出力する．
まずは，（なくてもよいが）後に出力されるデータ列の名前を出力し，
次に反復処理で，t(i),x(i),y(i)の値をn回出力する．
データの間にある「','」はデータ間の区切りを示すために入れている．「'　'」に挟まれた文字は，文字列として出力される．
いろいろなやり方があるが，ここではコンマ区切りとしている．

実行
==========

上記のプログラムコードをコンパイルし，実行するとターミナルに計算結果が表示される．

::

    t, x, y
  0.0000000000000000       ,   0.0000000000000000      ,   0.0000000000000000                                                                                                           
  0.10000000000000001      ,   0.86602540378443882     ,   0.49999999999999994                                                                                                           
  0.20000000000000001      ,   1.7320508075688776      ,   0.90199999999999991                                                                                                           
  0.30000000000000004      ,   2.5980762113533165      ,   1.2059999999999997                                                                                                           
  0.40000000000000002      ,   3.4641016151377553      ,   1.4119999999999997                                                                                                           
  0.50000000000000000      ,   4.3301270189221945      ,   1.5199999999999996                                                                                                           
  0.59999999999999998      ,   5.1961524227066338      ,   1.5299999999999994                                                                                                           
  0.69999999999999996      ,   6.0621778264910731      ,   1.4419999999999993                                                                                                           
  0.79999999999999993      ,   6.9282032302755123      ,   1.2559999999999991                                                                                                           
  0.89999999999999991      ,   7.7942286340599516      ,   0.97199999999999886   

コードの35~38行目に書いてある処理により時刻t，x座標，y座標の全ての計算結果が出力されている．
時刻tの列を見てみると分かりやすいが，0.1を足していくだけの計算であるはずなのに
「0.10000000000000001」のように小数点以下17桁目に値がある．
実数の計算精度に限界があることがわかる．

時間が経過するにつれて，x座標（2列目）は単調に増加，y座標（3列目）は一時増加の後，減少に転じているのが読み取れる．

.. index:: significant digits

捕捉（桁落ち）
====================

コンピュターが扱う数値には有効桁数に限界があることがわかった．
ただし注目している大きさに比べると非常に小さな誤差であり，あまり重要とは思わないかもしれない．
しかし，以下のような計算の場合には注意が必要である．

.. code-block:: fortran


    !+++++++++++++++++++++++++++++++++++++++++++++++
    program check_significant_digit
    !+++++++++++++++++++++++++++++++++++++++++++++++

    implicit none
    
    integer             i
    double precision    x1,x2,x3,y1,y2,y3
 
    x1 = 0.10000000000001d0
    x2 = 0.10000000000000d0
    x3 = x1-x2
    
    write(*,*) x1
    write(*,*) x2
    write(*,*) x3
    write(*,*) x3*1.0d14
    write(*,*) x1*1.0d14 - x2*1.0d14

    y1 = 1.0d10
    y2 = 1.0d-10
    y3 = y1+y2

    write(*,*) y1
    write(*,*) y2
    write(*,*) y3
    write(*,*) y3-y1

    end program check_significant_digit
    !+++++++++++++++++++++++++++++++++++++++++++++++

実行結果は以下の通り

::

    $ ./a.exe
    0.10000000000001000
    0.10000000000000001
    9.9920072216264089E-015
    0.99920072216264089
    1.0000000000000000

    10000000000.000000
    1.0000000000000000E-010
    10000000000.000000
    0.0000000000000000

``非常に近い数値の引き算(x1,x2)``

x2の小数点以下１７桁目に入った数値が，引き算によって誤差を多くしてしまい，
1.0d14倍したところでは，小数点以下4桁目がおかしくなってしまっている．
このような計算をしたいのならば，それぞれ1.0d14倍してから引き算した方がよい．

``大きさの極端に異なる数値(y1,y2)の足し算``

足したはずのy2がどこかにいってしまった．

（もう少し良い例はないか？）

.. warning::
    大きさの非常に近い数値の引き算や，大きさの極端に異なる数値の足し算が必要な場合は計算手順を工夫しよう．また誤差を軽減するためにも倍精度型実数を使うのがよい．