Pythonで差分法を実装して微分方程式を解いてみる

先日書籍を購入したので、その中で自分の為に残そうと思った記事を書いていきます。購入した書籍は以下です。

第一弾は微分方程式をパソコンに解かせるという話題です。運動方程式とか、波動方程式とか、マクスウェル方程式などの、身近な物理現象を記述してくれる式は全て微分方程式です。
大学生の時は、解が解析的に求まる状況ばかりを考えていましたが、現実が綺麗な系であることは稀です。
解が解析的に求まらなくても、数値的に微分方程式を解ける事がもあります。
今回の記事では、微分の差分化を使って簡単な運動方程式の解を得ます。

差分化とは
運動方程式
Python による実装
まとめ

差分化とは

差分化の話をする前に、1変数関数の微分の定義を思い出しましょう。関数$y =f(x) $の微分は以下の式で定義されます。
$$\begin{eqnarray}
\frac{dy}{dx} = \lim _{t \rightarrow 0} \frac{f(x+t) -f(x)}{t}
\end{eqnarray}$$
この式の意味は、$t $をどんどん小さくして行った時に、何か一つの値に収束していたら、それを$ \frac{dy}{dx} $と呼びます。という意味です。
世の中には、既に微分が出来ると分かっている関数が沢山存在します。つまり、このような関数では、常識で$t$を小さくしていけば、ある値に収束するという事です。小さい値$\Delta t $を与えるだけで大体の近似値を得る事が出来、誤差は$ (\Delta t) ^2 $程度です。
パソコンに極限を計算させることは出来ないかもしれませんが、人間が小さい数字$ \Delta t $ を与えて、以下の量を計算することは出来ます。
$$\begin{eqnarray}
\frac{dy}{dx} \sim \frac{f(x+\Delta t) -f(x)}{\Delta t}
\end{eqnarray}$$
この式を微分と思う事を微分の差分化と呼びます。

運動方程式

今回考える運動方程式は自由落下と、速度に比例した抵抗力がかかる場合の2種類です。運動方程式の説明は特にしませんが、例えば以下のサイトには学部で習うような物理の事が大体書いているのでお勧めです。具体的な計算がしたい人は別のサイトを見ないといけませんが。

https://eman-physics.net/

運動方程式は下記の2種類です。
$$\begin{eqnarray}
\frac{dy^2}{d^2t} &=&g \\
\frac{dy^2}{d^2t} &=& g – k \frac{dy}{dt}
\end{eqnarray}$$
質量は両辺割ったと思って書いてません。定数$g,k $は重力加速度と、速度に比例する抵抗力の比例係数です。どちらの式も解析的な解¹が求まります。
速度$ v= \frac{dy}{dt} $ については、以下のようになります。²
$$\begin{eqnarray}
v&=& – kt +a \\
v &=& \exp(-kt) + \frac{g}{k} +a
\end{eqnarray}$$
上式から、速度の大きさは抵抗力がある場合には $ v_{\infty} = \frac{g}{k} +a $ に収束することが分かります。
yについては以下のように書き表せます。
$$\begin{eqnarray}
y&=& – \frac{1}{2}t^2 +at +b\\
y &=& \frac{g}{k} t + a \exp(-kt) +b
\end{eqnarray}$$
運動方程式を差分化しましょう。まずは速度$ v= \frac{dy}{dt} $について差分化します。抵抗力のある場合だけ書いておいて、$ k=0 $とすれば自由落下の場合を再現出来るので、一つの場合だけ書きます。
$$\begin{eqnarray}
v(t+\Delta t) &=&v(t) + g\Delta t – kv(t) \Delta t
\end{eqnarray}$$
yについては、$ v = \frac{dy}{dt} $なので、いつも簡単な式になります。
$$\begin{eqnarray}
y(t+\Delta t) &=& y(t) + v(t) \Delta t
\end{eqnarray}$$
これで運動方程式をパソコンに解かせる準備が整いました。pythonでコードを書いて解いてみましょう。

Python による実装

人間が設定する必要のある事は、$ \Delta t $の値と、$ k$の値、そして初期値 $ y_0 , v_0 $ ³ です。コードは以下のようになると思います。
y軸を上側を正に取りたかったので、符号が変わってますが気にしないでください。


def parabolic(y, vel_y, stride, k):
""" 
y 位置
vei_y 速度
stride Δt
k 抵抗力の比例定数
"""
    g = 9.8
    y = y + vel_y * stride
    vel_y = vel_y - g * stride -k* vel_y *stride
    return [y, vel_y]

def parabolic_motion(y, vel_y, stride, k ):
    result =[] 
    cur = [y, vel_y]
    while(0<= cur[0]  and len(cur)<1000 ):#一生計算が続いても困るので条件をつけてます。
        result.append(cur)
        cur = parabolic(cur[0], cur[1], stride, k)
    return result

$ (y_0, v_0 )=(1000, 10) , k=0.1, \Delta t =0.3 $, として解かせると、以下のようなグラフになります。

パット見た感じ良さそうです。自由落下の場合は速度がどんどん大きくなっていくのに対して、抵抗力がある場合はある速度に収束していく様子が再現出来ました。大体-80から-90くらいに収束しているように見えますが、今回与えた定数では$ v_{\infty} =- \frac{g}{k} +a =88 $となっているので、大体理論値と一致しています。
別の記事で、偏微分方程式の場合も解いてみたいと思います。