Python で ode45()関数を模倣する

MATLAB の ode45() 関数は、適応型ルンゲ＝クッタ法で常微分方程式（ODE）の初期値問題を解きます。Python には文字どおり ode45() という名前の関数はありませんが、SciPy の solve_ivp() が同じワークフローをカバーします。方程式を記述し、ソルバーに時間区間と初期状態を渡し、返されたオブジェクトから解を読み取ります。

MATLAB の ode45() には SciPy の solve_ivp() を使う

scipy.integrate の solve_ivp() 関数は、ode45() に最も近い実用的な代替です。デフォルト手法の RK45 は、自動刻み幅制御を備えた次数 5(4) の陽的ルンゲ＝クッタ公式であり、MATLAB で ode45() が使うのと同じアルゴリズム系です。そのため、解の変化が緩やかな区間では大きなステップを、変化が速い区間では小さなステップを取ります。

ただし、そのまま置き換えられるわけではありません。引数名が異なり、結果も t と y の配列のペアではなく 1 つのオブジェクトとして返されます。それでも、導関数・時間区間・初期状態という全体のパターンは両者の間で直接対応します。SciPy がまだ利用できない場合は、pip install scipy でインストールしてください。

ODE を定義して solve_ivp() を呼び出す

基本的な呼び出しは solve_ivp(fun, t_span, y0, method="RK45", t_eval=None, args=None, ...) です。主な入力は次のとおりです。

• fun: 現在の時刻 t と状態 y に対する導関数 dy/dt を返す呼び出し可能オブジェクト fun(t, y)。単一の方程式では 1 つの値を返し、連立方程式では状態成分ごとの導関数を含む配列を返します。
• t_span: 積分区間の開始と終了を表す (t_start, t_end) のペア。
• y0: t_start における状態。その長さが方程式系の本数を決めます。
• method（省略可）: 積分スキーム。デフォルトは 'RK45' です。
• args（省略可）: fun(t, y, *args) のように関数へ渡される追加の値。
• rtol と atol（省略可）: 適応刻み幅を制御する相対許容誤差と絶対許容誤差。

次の例は、初期値 y(0) = 1 から区間 [0, 2] で dy/dt = 2y を解きます。この呼び出しは結果オブジェクトを返し、その t 属性と y 属性を表示すると、ソルバーが選んだ時点と各時点での状態がわかります。次のコードと出力に表示されるラベルや文字列は英語のままです:

from scipy.integrate import solve_ivp


def growth(t, y):
    return 2 * y


result = solve_ivp(growth, (0, 2), [1])
print("t:", result.t)
print("y:", result.y[0])

t: [0.         0.07580098 0.56590702 1.18757739 1.85907272 2.        ]
y: [ 1.          1.16369695  3.10133569 10.75272756 41.18553629 54.59497095]

報告される時刻の間隔が不均一なのは、solve_ivp() が一定の増分で進むのではなく、各ステップを解の局所的な振る舞いに適応させるためです。

ソルバーの結果を読み取る

solve_ivp() は scipy.integrate.OdeResult オブジェクトを返します。最もよく使う属性は次のとおりです。

• t: 報告された時点を格納する NumPy 配列。t_eval を指定しない場合、これらはソルバー内部の適応ステップです。
• y: 形状 (n, len(t)) の 2 次元 NumPy 配列。行 i には報告された各時刻での状態の成分 i が入るため、y[:, k] は t[k] における状態です。単一の方程式では、y[0] が解曲線です。
• sol: result.sol(t) として呼び出し、区間内の任意の時刻で解を評価できる連続的な OdeSolution 補間オブジェクト。dense_output=True を渡した場合にのみ設定され、それ以外では None です。
• success と message: 積分が完了したかどうかを示すブール値と、どのように終了したかを説明する人間が読めるメッセージ。
• status: 区間が完了した場合は 0、終了イベントが発生した場合は 1、ソルバーが進展できなかった場合は -1。

失敗した実行でも部分的な配列が返されることがあるため、y を使う前に success（または status）を確認する習慣を付けるとよいでしょう。

t_eval で出力時刻を指定する

適応ステップの時刻は、報告したい時刻と一致するとは限りません。ソート済みの時刻配列を t_eval に渡すと、solve_ivp() は内部ステップではなく、指定した時刻ちょうどの解を返します。次のコードと出力に表示されるラベルや文字列は英語のままです:

from scipy.integrate import solve_ivp


def growth(t, y):
    return 2 * y


times = [0.5, 1.0, 1.5, 2.0]
result = solve_ivp(growth, (0, 2), [1], t_eval=times)
print("t:", result.t)
print("y:", result.y[0])

t: [0.5 1.  1.5 2. ]
y: [ 2.71823524  7.38683161 20.07910026 54.59497095]

これで t には要求した時刻だけが含まれ、y にはそれぞれの時刻に補間された状態が格納されます。t_eval の値はすべて t_span の内側になければなりません。積分区間の外にある時刻を要求すると、その値は t_span 内になければならないというエラーで呼び出しが停止します。端点は許可されているため、t_eval には t_start と t_end の両方を含められます。