pythonでルンゲクッタ法を実装した.
レスラー方程式を与えて描画する.
プログラム内関数のdxdt, dydt, dzdtを書き換えれば別の微分方程式も解ける.
下記プログラムはローレンツ方程式を解くようにも書き換えられる.
import numpy as np def dxdt(a, b, c, x, y, z): # 引数は,「関数のパラメータを並べたもの」,「変数を並べたもの」の順番で並べる.他の連立微分方程式と同じ形式で引数を並べる. return -y - z def dydt(a, b, c, x, y, z): return x + a * y def dzdt(a, b, c, x, y, z): return b + x * z - c * z def _handler(func, *args): return func(*args) pdelist = [dxdt, dydt, dzdt] # 連立微分方程式を追加する場合は,関数名をこのlistにも追加すること. def runge_kutta(funcs, param, vold, delta): """ ルンゲクッタ法で時間を含まない微分方程式を解く :param funcs: 関数のlist.上記のpdelistのように並べる. :param param: パラメータを並べたndarray(np.array) :param vold: 現在時刻における未知関数の値を並べたndarray.pdelistに含まれる関数の引数のx, y, zに相当する. :param delta: 計算の1ステップあたりに進める時間 :return: 計算を1ステップ進めた後の時刻における未知関数の値を並べたndarray.pdelistに含まれる関数の引数のx, y, zに相当する. """ k1 = [] k2 = [] k3 = [] k4 = [] for fnc in funcs: k1.append(_handler(fnc, *param, *vold)) k11 = vold + delta * 0.5 * np.array(k1) for fnc in funcs: k2.append(_handler(fnc, *param, *k11)) k22 = vold + delta * 0.5 * np.array(k2) for fnc in funcs: k3.append(_handler(fnc, *param, *k22)) k33 = vold + delta * np.array(k3) for fnc in funcs: k4.append(_handler(fnc, *param, *k33)) return vold + (np.array(k1) + 2 * np.array(k2) + 2 * np.array(k3) + np.array(k4)) * delta / 6 def main(): import pickle import matplotlib.pyplot as plt param = np.array([0.2, 0.2, 7.1]) # レスラー方程式のパラメータ init = np.array([4, 4, 1]) # 微分方程式の初期値 delta = 0.1 # 1ステップあたりの経過時間 result = [] vnew = init for i in range(1, 10000): result.append(vnew) vnew = runge_kutta(pdelist, param, vnew, delta) result.append(vnew) roessler = np.array(result) # 解の変数.roessler[時刻(sample), 変数の番号] with open("signal_roessler.pickle", mode="wb") as f: # 解の変数を保存 pickle.dump(roessler, f) plt.plot(roessler[:, 0], roessler[:, 1]) # 解の変数の0番目(レスラー方程式のx)を横軸に,解の変数の1番目(レスラー方程式のy)を縦軸に取り,プロットする. plt.xlabel("x = roessler[:, 0]") plt.ylabel("y = roessler[:, 1]") if __name__ == "__main__": import matplotlib.pyplot as plt main() plt.show()
上記プログラムを実行すると以下のグラフが表示される.