4.1. 常微分方程式の解法

4.1.1. 概観

常微分方程式（Ordinary Differential Equation; ODE）を考えます。 一変数の最も簡単な方程式は以下の形で与えられます：

\frac{dy(t)}{dt} = f(t, y(t))

2階以上の微分を含む場合は、新しい変数を導入することで連立1階微分方程式に書き換えられます。 例えば、ニュートン方程式 m d^2x(t)/dt^2=F(t) の場合、速度 v(t) を導入すると dx(t)/dt=v(t) と m dv(t)/dt=F(t) の2つの1階微分方程式に分けられます（詳細は 【例題】ニュートン方程式 を参照）。 そこで、\bm{y}=(y_1, y_2, \cdots, y_N) と定義すると、ODEの一般的な表式は

(1)\frac{d\bm{y}(t)}{dt} = \bm{f}(t, \bm{y}(t))

と表せます。ここで、太字の \bm{f} はN個の関数をベクトルとしてまとめたものを表します。 この方程式を初期条件 \bm{y}(0)=\bm{y}_0 の下で解く問題を考えます。

ODEを数値的に解く場合には変数 t を離散化し、t=t_i から t=t_{i+1} へ \bm{y}(t) を逐次的に更新していきます：

(2)\bm{y}(t_{i+1}) = \bm{y}(t_i) + \Delta \bm{y}

\Delta \bm{y} を見積もる方法として以下のようなものがあります：

  • 陽解法 (Explicit method)

  • 陰解法 (Implicit method)

  • 予測子・修正子法 (Predictor-Corrector method)

陽解法 は最も直接的な方法で、その代表が Euler法 です。 (1) 式の微分を1次の前進差分で置き換えて、\Delta \bm{y}= f(t_i, \bm{y}(t_i)) \Delta t とします。 これは (\Delta t)^2 を無視する近似ですので、精度は \mathcal{O}(\Delta t) です。 実用性はありません。 精度と実装の簡単さのバランスが取れていて最もよく使われるのが Runge-Kutta法 です。 公式はいくつかありますが、\mathcal{O}((\Delta t)^4) の公式が最も有名です。 単にRunge-Kutta法という場合には4次の公式を指します。

陰解法は (1) 式の微分を後退差分で置き換えたものです。 1次の公式の場合、 \Delta \bm{y}= f(t_{i+1}, \bm{y}(t_{i+1})) \Delta t となります。 \bm{y}(t_{i+1}) を求めるために \bm{y}(t_{i+1}) を使うというのは一見不思議なようですが、 (2) 式に \Delta \bm{y} を代入して得られる方程式を \bm{y}(t_i) について解くことで計算できます。 陽解法よりも安定ですが、各ステップでこの方程式を解くため、陽解法よりもだいぶ複雑になります。 この解説では取り上げません。

予測子・修正子法 は陽解法と陰解法を組み合わせた方法です。 まず陽解法で \bm{y}(t_{i+1}) を求めておき（これを予測子と呼ぶ）、陰解法を使って予測子を修正します。 予測子を用いることで、陰解法の方程式を短い時間で解けるようにするという見方もできます。 ただ、[Press] によると

予測子・修正子法はもう時代遅れで、常微分方程式のほとんどの問題について、もはや最良の方法ではないのではなかろうか。 高精度を要したり右辺の計算が厄介な場合は、Bulirsch-Stoer方が良い。 使いやすさでは、あるいは低精度でよいなら、適応刻み幅のRunge-Kutta法に軍配が上がる。

とあります。 この解説ではライブラリを使う方針ですので参考までに。

4.1.2. 不安定性

ODEの数値解法では、刻み幅 \Delta t をある程度小さく取らないと、計算が不安定になる（解が発散したりする）場合があります。 この不安定性はアルゴリズム特有のもので、バグではありません。 もし、数値計算でおかしな結果が得られたら、\Delta t を十分小さく取って、解が変わるかを確認してください。 これで解消されればバグではありません。 実際に \Delta t をどの程度まで小さくすればよいかは、解法や方程式に依存します。

4.1.3. 適応刻み幅制御

不安定性を回避するため、また、精度を上げるために、刻み幅 \Delta t はなるべく小さく取りたい。 しかし \Delta t を小さくすると多くのステップ数が必要になり、計算時間がかかってしまう。 このようなジレンマがあります。

適度な \Delta t の大きさは、一般に、計算してみないと分かりません。 また、時刻 t にも依存します。 \bm{y}(t) の変化が激しければ \Delta t を小さくする必要がありますが、変化が緩やかであれば \Delta t は大きくできます。

そこで、\Delta t を自動的に調整する方法が考案されています。 これを、適応刻み幅制御 （adaptive stepsize control）と言います。 ODEに依らず、関数の説明で adaptive と書いてあったら、連続変数の離散化などが自動的に最適化されることを意味しています。 adaptiveな手法の場合は、刻み幅 \Delta t ではなく、精度を指定します。

注釈

adaptiveな手法では、内部で精度を見積もり、精度が足りていなければ \Delta t を小さく、十分足りていれば \Delta t を大きく変えながらステップを進めます。精度の見積もりかたは実装依存ですが、例えば、\Delta t ×1ステップの計算と \Delta t/2 ×2ステップの計算を比較して、その差から精度を見積もります。 例えば、[Press] に解説があります。

4.1.4. scipy.integrate.solve_ivpの使い方

常微分方程式の数値解法はSciPyの scipy.integrate モジュールに含まれています。 scipy.integrate.solve_ivp 関数を使います。

注釈

他にも古い関数（Old API）として、 関数版の scipy.integrate.odeint と クラス版の scipy.integrate.ode があります。 インターフェース（関数の使い方）が古い作りになっていますが、非推奨というわけではないようです。

While the interface to them is not particularly convenient and certain features are missing compared to the new API, the solvers themselves are of good quality and work fast as compiled Fortran code. In some cases, it might be worth using this old API.

この2つであれば、ode の方が機能が豊富です。

さて、scipy.integrate.solve_ivp の使い方を見ていきます。 関数の引数は以下のようになっています：

scipy.integrate.solve_ivp(fun, t_span, y0, method='RK45', t_eval=None,
                          dense_output=False, events=None, vectorized=False,
                          args=None, **options)

fun は方程式 (1) の f(t, \bm{y}) に対応する関数です。 fun(t, y) のように定義します。t は実数、y はNumPy配列です。 fun(t, y, a, b) のように、他に引数を追加しても構いません。方程式に含まれるパラメーターの指定などに使います。 a, b の値は args 引数を使って args=(a0, b0) のように与えます。 t_span は時刻の最初と最後を t_span=(0, 10.0) や t_span=[0, 10.0] のようにタプル形式またはリスト形式で与えます。y0 は初期条件です。

method でアルゴリズムを指定します。 以下のアルゴリズムが選択可能です（2020年10月時点）。

  • ‘RK45’ (default): Explicit Runge-Kutta method of order 5(4).

  • ‘RK23’: Explicit Runge-Kutta method of order 3(2).

  • ‘DOP853’: Explicit Runge-Kutta method of order 8.

  • ‘Radau’: Implicit Runge-Kutta method of the Radau IIA family of order 5.

  • ‘BDF’: Implicit multi-step variable-order (1 to 5) method based on a backward differentiation formula for the derivative approximation.

  • ‘LSODA’: Adams/BDF method with automatic stiffness detection and switching.

上記の全てにおいて、adaptive stepsize controlが実装されています。 通常はデフォルトのまま（'RK45'）で問題ないと思います。 詳しい説明は公式ドキュメント scipy.integrate.solve_ivp を見てください。

最後の **options は method に依存したオプションを渡すことを意味します。 例えば、RK45で利用可能なオプションは scipy.integrate.RK45 クラスのドキュメントを見ると分かります。

class scipy.integrate.RK45(fun, t0, y0, t_bound, max_step=inf,
                           rtol=0.001, atol=1e-06,
                           vectorized=False, first_step=None, **extraneous)

max_step 以降が指定可能です。 この中で重要なものは rtol と atol で、それぞれ相対誤差、絶対誤差を表します。 必要に応じて、精度を上げたい場合は小さく、時間を節約したい場合は大きくします。