オイラー法でボール投げ

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

運動方程式

運動について何か知りたければ，まずはとにかく運動方程式である．

水平方向座標を \(x\) ，高さ方向座標を \(y\) とし，時間 \(t\) とする．ボールの質量が \(m\) で重力加速度を \(g\) のとき，

運動方程式は以下のように書ける．

(1)\[\begin{split}m\frac{d^2 x}{dt^2} &= 0 \\ m\frac{d^2 y}{dt^2} &= -mg\end{split}\]

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

(2)\[\begin{split}x &= x_0 + v_0\cos{\theta}\ t \\ y &= y_0 + v_0\sin{\theta}\ t - \frac{1}{2}gt^2\end{split}\]

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

数値積分（オイラー法）

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

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

(3)\[\dot f = \Phi (t)\]

左辺の微分演算は，

\[\dot f = \frac{df}{dt} = \lim\limits_{\epsilon \to 0} \frac{f(t+\epsilon)-f(t)}{\epsilon}\]

であるが，これを

\[\lim\limits_{\epsilon \to 0} \frac{f(t+\epsilon)-f(t)}{\epsilon} \approx \frac{f(t+\Delta t)-f(t)}{\Delta t}\]

とする． \(\Delta t\) が小さいのであれば，このように近似してもよいだろう．

すると，元の式は，

\[\dot f = \frac{f(t+\Delta t)-f(t)}{\Delta t} = \Phi (t)\]

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

(4)\[f(t+\Delta t) = f(t) + \Phi (t) \Delta t\]

この式は， \(f(t)\), \(\Phi (t)\), \(\Delta t\) が分かっていれば， \(f(t+\Delta t)\) が計算できることを意味する．

つまり，時刻 \(t\) での情報 \(f(t)\) , \(\Phi (t)\) と時間刻み \(\Delta t\) により，時刻 \(t + \Delta t\) での情報 \(f(t + \Delta t)\) が得られる．

この計算方法をオイラー陽解法（Euler’s explicit scheme）という．

Note

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

では，計算してみよう

計算手順

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

Note

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

すると，式 (1) の二つの式は以下のように書ける．

\[\begin{split}\dot x &= u, \dot u &= 0 \\ \dot y &= v, \dot v &= -g\end{split}\]

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

(5)\[\begin{split}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\end{split}\]

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

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

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

Fortranで計算しよう

コンピュータは単純な計算を繰り返し行うことが得意である．

Fortranを使って計算手順を記述し，コンピュータで計算を実行してみよう．

!+++++++++++++++++++++++++++++++++++++++++++++++
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になれていると，面倒くさいと思ってしまうが仕方がない．

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

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

program ball_throw


end program ball_throw

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

implicit none

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

integer                     i
integer,parameter::         n=10

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

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⁰ =9.8である．gは重量加速度であり，変化する値ではないのでparameterがついている． 次のpi = acos(-1.0d0)では，piという変数にacos(-1.0d0)で計算した値を入れている．このように計算した結果を使うこともできる． acos(x)はcos^-1 (x)を計算する関数である． つまりpiにcos^-1 (-1)の計算結果を入れている．円周率 \(\pi\) を計算して求めている． もちろんpiに値をいれるのにpi=3.141592....d0と数値を使ってもよい．

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

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

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をラジアン単位に変換している．

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

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

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回の反復計算を行う． 反復計算の中身は．．．じっと見ていればわかると思うが， (5) の計算である． 「=」の右側（例えば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)を計算する．

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列目）は一時増加の後，減少に転じているのが読み取れる．

捕捉（桁落ち）

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

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

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