今回は、Pythonによる常微分方程式を解く際によく使われるscipyライブラリのメモである。
Scipyによる常微分方程式
- Pythonを用いて常微分方程式を解く際は、ライブラリを用いるのが最も手軽で、かつ高精度な結果が得られる。
- 常微分方程式を解くライブラリとして、Scipyがある。
- Scipyの中のscipy.integrateモジュールに含まれているslove_ivp関数を使う。solve_ivpのivpは"Initial Value Problem"の略である。
- Pythonコード中では、以下のように呼び出す。
from scipy.integrate import solve_ivp
- 使い方は、以下のようなコーディングの形となる。
scipy.integrate.solve_ivp(fun, t_span, y0, methid='RK5', t_eval=None, dense_output=False, events=None, vectorized=False, arags=None, **options)
solve_ivp関数の引数は次の表に示す。
| 引数 | 説明 |
|---|---|
| fun(関数オブジェクト) | 微分方程式の右辺$\boldsymbol{f}(t,\boldsymbol{y})$に対応する関数。$fun(t,y)$のように定義した関数を与える。$fun(t,y,a,b)$のように、他に引数を追加ししても構わない。方程式に含まれるパラメータの指定などに使う。この場合、$a$、$b$の値はargs引数を使って、args=(1.0,-3.5)のように与える。 |
| t_span(tuple) | 時刻の最初と最後をt_span=(0,10.0)やt_span=[0,10.0]のようにタプルやリストの形式で与える。 |
| y0(ndarray(n,)) | 初期値をNumpy配列として与える。 |
| method(str) | アルゴリズムを指定する。 |
| dense_output(bool) | 任意の時刻における結果を得るには、Trueにしておく必要がある。 |
| events(関数オブジェクト) | ある条件を満たす時刻tを求めたい場合に使用する。例えば、球が地面に到達する時刻を知りたい場合など。 |
| args(tuple) | funに与えた関数がパラメータ(3つ目以降の引数)をとる場合に、パラメータの値を入力する。 |
| 戻り値 | 説明 |
| t(np.ndarray(n_points,)) | 自動的に選ばれた時刻tの配列 |
| y(np.ndarray(n,n_points)) | 時刻tにおける結果$\boldsymbol{y}(t)$の配列。 |
| sol(関数オブジェクト) | sol(t_dense)のように呼び出すことにより、任意の時刻t_denseにおける結果が得られる(dense_ouput=Trueの場合のみ) |
| message(str) | 計算終了時の状態を表す短い文 |
次にmethodに指定できるアルゴリズムを次の表に示す。なお、アルゴリズムの説明は以下のscipyの公式マニュアルを参照している。
https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html
| オプション | アルゴリズム |
|---|---|
| 'RK5' | (default)Explicit Runge-Kuttta method of order 5(4) |
| 'RK2' | 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 diffenciation formula for the derivative approximation |
| 'LSODA' | Adams/BDF method with automatic stiffness detection and switching |
