Ccmmutty logo
Commutty IT
・0 pv・9 min read

微分方程式の数値解析の基本①

数学数値解析Python
https://cdn.magicode.io/media/notebox/blob_jh4ygJd
昔ちょっとやったことがある微分方程式の数値解析をまた1から取り組んでみたくなったので復習を兼ねて色々やってみる。
記事タイトルに①とついているが続くかは自分のやる気次第。
少しだけ複雑な方程式でやってみようかと思ったが参考資料通りにやろうとしても結構苦戦したので、まず簡単な自由落下を対象にしながら整理していくことにした。
数値解析の基本中の基本とも言われるオイラー法と4次のルンゲ・クッタ法の2つで解いてみて両者の結果を比較する。

自由落下の運動方程式と理論解

おそらく最も単純な物理現象の一つである自由落下の運動方程式は、yy座標を鉛直上向を正に取り重力加速度をg[m/s2]g[m/s^2]、落下する物体の質量をm[kg]m[kg]とすれば
md2ydt2=−mgm \dfrac {d^2y}{dt^2} = -mg
となるので、解析対象の微分方程式は
d2ydt2=−g\dfrac{d^2y}{dt^2} = -g
この微分方程式の両辺を時間ttで積分すれば速度の関数v(t)v(t)が得られる。（ただしv(0)=v0v(0)=v_0）
dydt=−gt+v0=v(t)\dfrac{dy}{dt} = -gt + v_0 = v(t)
さらに両辺をttで積分すれば位置の関数y(t)y(t)が得られる。（ただしy(0)=y0y(0)=y_0）
y=−12gt2+v0t+y0=y(t)y = -\dfrac{1}{2}gt^2 + v_0t + y_0 = y(t)
このように自由落下の微分方程式は手計算で簡単に速度や位置の関数を求めることができるのでわざわざ数値的に解くまでもないが勉強のためにやる。

オイラー法とルンゲ・クッタ法について

オイラー法

xxの関数yyをxxで微分した1階の常微分方程式の一般形は下記のように表される。
dydx=f(x,y)\dfrac{dy}{dx} = f(x, y)
ここで、hhを微小な刻み幅とするとy(x+h)y(x+h)の値はテイラー展開で下記のように表せる。
y(x+h)=y(x)+dydxh+12d2ydx2h2+…y(x + h) = y(x) + \dfrac{dy}{dx}h + \dfrac{1}{2}\dfrac{d^2y}{dx^2}h^2 + …
特に、2階以上の微分項をまとめてO(h2)O(h^2)と書くと
y(x+h)=y(x)+dydxh+O(h2)y(x + h) = y(x) + \dfrac{dy}{dx}h + O(h^2)
オイラー法はこれの1階微分の項までを使いy(x)y(x)の値から次の計算ステップの値y(x+h)y(x+h)を計算する。
y(x+h)=yi+1, y(x)=yi, dydx∣x=xi=f(xi,yi)y(x + h) = y_{i+1},\space y(x) = y_i,\space \left. \dfrac{dy}{dx}\right| _{x=x_i} = f(x_i, y_i)とすれば
yi+1=yi+f(xi,yi)hy_{i+1} = y_i + f(x_i, y_i)h
のような漸化式で次ステップのyi+1y_{i+1}を計算できる。
ただし、この漸化式にはO(h2)O(h^2)の部分が含まれていないので、計算が繰り返されればその分だけ誤差が蓄積していくことになる。
O(h2)O(h^2)を無視することによる誤差を打ち切り誤差といい、hhが大きすぎると打ち切り誤差の影響が大きくなり精度が下がってしまい、かといってhhを小さくしすぎると計算時間が増加してしまうことになる。

4次のルンゲ・クッタ法

hhを小さめに取らなくても精度良く計算する方法の1つにルンゲ・クッタ法がある。
ルンゲ・クッタ法と言ってもその種類はたくさんあり、中でも4次のルンゲ・クッタ法という方法が良く利用されるらしい。
この方法では下記のように次々とk1k_1〜k4k_4を計算してyi+1y_{i+1}を求める。
k1=f(xi,yi)hk2=f(xi+h2,yi+k12)hk3=f(xi+h2,yi+k22)hk4=f(xi+h,yi+k3)hyi+1=yi+k1+2k2+2k3+k46\begin{aligned} k_1&=f(x_i, y_i)h \\[0.5em] k_2&=f(x_i + \dfrac{h}{2}, y_i + \dfrac{k_1}{2} )h \\[0.5em] k_3&=f(x_i + \dfrac{h}{2}, y_i + \dfrac{k_2}{2} )h \\[0.5em] k_4&=f(x_i + h, y_i + k_3)h \\[0.5em] y_{i+1}&= y_i + \dfrac{k_1 + 2k_2 + 2k_3 + k_4}{6} \end{aligned}
k1k_1はオイラー法におけるyyの増分と同じである。
なお、参考書や参考記事によっては
k1=f(xi,yi)k2=f(xi+h2,yi+k1h2)k3=f(xi+h2,yi+k2h2)k4=f(xi+h,yi+k3h)yi+1=yi+k1+2k2+2k3+k46h\begin{aligned} k_1&=f(x_i, y_i) \\[0.5em] k_2&=f(x_i + \dfrac{h}{2}, y_i + \dfrac{k_1h}{2} ) \\[0.5em] k_3&=f(x_i + \dfrac{h}{2}, y_i + \dfrac{k_2h}{2} ) \\[0.5em] k_4&=f(x_i + h, y_i + k_3h) \\[0.5em] y_{i+1}&= y_i + \dfrac{k_1 + 2k_2 + 2k_3 + k_4}{6}h \end{aligned}
のようにhhを付ける位置が若干異なっている記載もあるが、knk_nをyyの増分と見るか（前者）yyの勾配と見るか（後者）の違いなだけで行き着く先は同じである。
この記事では前者の書き方で統一している。

自由落下の運動に適用する

オイラー法、4次のルンゲ・クッタ法は1階の常微分方程式を解くための方法だが、今計算したい自由落下の微分方程式は
d2ydt2=−g\dfrac{d^2y}{dt^2} = -g
という2階の常微分方程式である。
なので、媒介変数vvを使って2階の常微分方程式を1階の連立微分方程式に分けてやれば上記の方法を適用して数値的に解くことができる。
dydt=v=fy(t,y,v)dvdt=−g=fv(t,y,v)\begin{aligned} \dfrac{dy}{dt}&= v = f_y(t, y, v)\\[0.5em] \dfrac{dv}{dt}&= -g = f_v(t, y, v) \end{aligned}
このようにすればオイラー法の漸化式は下記のようになる。
hhは計算ステップの間隔、つまり計算の時間刻み幅になる。
vi+1=vi+dvdth=vi−gh=vi+fv(ti,yi,vi)hyi+1=yi+dydth=yi+vh=yi+fy(ti,yi,vi)hv_{i+1} = v_i + \dfrac{dv}{dt}h =v_i - gh =v_i + f_v(t_i, y_i, v_i)h \\[0.5em] y_{i+1} = y_i + \dfrac{dy}{dt}h =y_i + vh =y_i + f_y(t_i, y_i, v_i)h
4次のルンゲ・クッタ法の計算も同様で、速度を計算するためのk1k_1〜k4k_4と位置を計算するためのl1l_1〜l4l_4を順番に計算する。
k1=fv(ti,yi,vi)hl1=fy(ti,yi,vi)hk2=fv(ti+h2,yi+l12,vi+k12)hl2=fy(ti+h2,yi+l12,vi+k12)hk3=fv(ti+h2,yi+l22,vi+k22)hl3=fy(ti+h2,yi+l22,vi+k22)hk4=fv(ti+h,yi+l3,vi+k3)hl4=fy(ti+h,yi+l3,vi+k3)hvi+1=vi+k1+2k2+2k3+k46yi+1=yi+l1+2l2+2l3+l46\begin{aligned} k_1&=f_v(t_i, y_i, v_i)h \\[0.5em] l_1&=f_y(t_i, y_i, v_i)h \\[1em] k_2&=f_v(t_i + \dfrac{h}{2}, y_i + \dfrac{l_1}{2}, v_i + \dfrac{k_1}{2})h \\[0.5em] l_2&=f_y(t_i + \dfrac{h}{2}, y_i + \dfrac{l_1}{2}, v_i + \dfrac{k_1}{2})h \\[1em] k_3&=f_v(t_i + \dfrac{h}{2}, y_i + \dfrac{l_2}{2}, v_i + \dfrac{k_2}{2})h \\[0.5em] l_3&=f_y(t_i + \dfrac{h}{2}, y_i + \dfrac{l_2}{2}, v_i + \dfrac{k_2}{2})h \\[1em] k_4&=f_v(t_i + h, y_i + l_3, v_i + k_3)h \\[0.5em] l_4&=f_y(t_i + h, y_i + l_3, v_i + k_3)h \\[1em] v_{i+1}&= v_i + \dfrac{k_1 + 2k_2 + 2k_3 + k_4}{6} \\[0.5em] y_{i+1}&= y_i + \dfrac{l_1 + 2l_2 + 2l_3 + l_4}{6} \end{aligned}
なお、今回の自由落下の場合であればfvf_vは定数項−g-gしかないしfyf_yはvvをそのまま返すだけなのでわざわざ引数にt,yt, yあるいはvvを渡さなくても済むはずである。
しかし、もっと複雑な常微分方程式を対象にする場合はfvf_vやfyf_yにt,yt, yやvvの項が現れることもあるので、今回はあえて直接影響がない変数も引数として渡すような記述にしている。
この辺りは参考にした資料によっては使われない引数が省略されていたりされていなかったり記述に差がありちょっとややこしいところだったので理解に時間がかかった。

数値計算してみる

外部入力で物体の鉛直方向の初速度と初期位置を入力し、オイラー法による数値解、4次のルンゲ・クッタ法による数値解、理論解でそれぞれ計算した結果をプロットして比べる。
時間刻み幅はh=0.01h=0.01に設定している。
速度と位置の両方を計算しているが、とりあえずプロットは位置だけにしている。
Pythonで計算プログラムを作るのであれば簡単に常微分方程式の数値解を得られるSciPyのodeintなどがあるのだがそれはまた別の機会に。
import matplotlib.pyplot as plt

G = 9.80665  # 重力加速度
h = 0.01  # 時間刻み幅


def f_v(t, y, v):
    """dv/dt"""
    return -G


def f_y(t, y, v):
    """dy/dt"""
    return v


def calc_euler_method(t, y, v):
    """オイラー法により次ステップのv, yを計算する"""
    v_next = v + f_v(t, y, v) * h
    y_next = y + f_y(t, y, v) * h

    return v_next, y_next


def calc_runge_kutta_method(t, y, v):
    """4次のルンゲ・クッタ法により次ステップのv, yを計算する
        速度vを計算する時の係数をk1〜k4、位置yを計算する時の係数をl1〜l4とする
    """
    k1 = f_v(t, y, v) * h
    l1 = f_y(t, y, v) * h
    k2 = f_v(t + h / 2, y + l1 / 2, v + k1 / 2) * h
    l2 = f_y(t + h / 2, y + l1 / 2, v + k1 / 2) * h
    k3 = f_v(t + h / 2, y + l2 / 2, v + k2 / 2) * h
    l3 = f_y(t + h / 2, y + l2 / 2, v + k2 / 2) * h
    k4 = f_v(t + h, y + l3, v + k3) * h
    l4 = f_y(t + h, y + l3, v + k3) * h
    v_next = v + (k1 + 2 * k2 + 2 * k3 + k4) / 6
    y_next = y + (l1 + 2 * l2 + 2 * l3 + l4) / 6

    return v_next, y_next


def calc_theoretical_solution(t, y0, v0):
    """理論解を計算する"""
    v = -G * t + v0
    y = -G * t ** 2 / 2 + v0 * t + y0

    return v, y


def freefall(y0, v0):
    """自由落下の計算"""
    t = 0.0  # 時刻
    t_list = [t]
    ve = vr = vt = v0
    ye = yr = yt = y0
    ye_list = [ye]  # yの数値解（オイラー法）
    yr_list = [yr]  # yの数値解（ルンゲ・クッタ法）
    yt_list = [yt]  # yの理論解
    while ye > 0 and yr > 0 and yt > 0:
        t += h
        ve, ye = calc_euler_method(t, ye, ve)
        vr, yr = calc_runge_kutta_method(t, yr, vr)
        vt, yt = calc_theoretical_solution(t, y0, v0)

        t_list.append(t)
        ye_list.append(ye)
        yr_list.append(yr)
        yt_list.append(yt)

    return t_list, ye_list, yr_list, yt_list


if __name__ == '__main__':
    y0 = float(input('初期高度y0:'))
    v0 = float(input('初速度v0:'))

    # 各種方法で自由落下を計算
    t_list, ye_list, yr_list, yt_list = freefall(y0, v0)

    # 位置のプロット
    plt.title('Freefall  (time interval: {}s)'.format(h))
    plt.plot(t_list, yt_list, label='theoretical',
             color='green', linestyle='None', marker='o', markersize=1.5)
    plt.plot(t_list, ye_list, label='Euler', color='blue', linewidth=1)
    plt.plot(t_list, yr_list, label='Runge-Kutta', color='red', linewidth=1)
    plt.xlabel("time [s]")
    plt.ylabel("y [m]")
    plt.legend()
    plt.show()
初速度50m/s50m/s（鉛直投げ上げ）、初期位置100m100mで計算してみる。
h=0.01h=0.01であればオイラー法、4次のルンゲ・クッタ法ともに理論解とほとんど一致しているが、時間刻み幅を10倍に拡大、つまりh=0.1h=0.1にするとオイラー法での計算の誤差が目立ってくる。

Discussion

コメントにはログインが必要です。