pythonでレスラー方程式の解軌道を描画（ルンゲクッタ法）

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()