scipy.integrate.solve_ivpについて浅く掘る
scipy.integrate.solve_ivpの実装について簡単に見ていきます。 Solving Ordinary Differential Equations I [1]と実際のコード[2]を参考にしています。
設定
今回はLorenz方程式を使用します。
Deterministic Nonperiodic Flow
$ \begin{cases} \dot{x}=\sigma(y-x)\\ \dot{y}=rx-y-xz\\ \dot{x}=xy-bz \end{cases} $
パラメータは代表的な$\sigma=10,b=8/3,r=28$を使います。
=
= - * + *
= - * + * -
= * - *
return
, ,
今回は初期値を$(0.1, 0.1, 0.1)$とし、時間$t=0$から$t=40$で解いていきます。
使うライブラリは以下の通りです。
solve_ivp
scipy.integrate.solve_ivp
を使って解いていきます。
方程式を解く際にscipy.integrate.odeint
やscipy.integrate.ode
を使っている例が多いですが、現在はOld API
に指定されています。
Integration and ODEs | Old API
対象の方程式、時間のスパン、初期値を与え、解きます。
=
計算結果はOdeResult
というオブジェクトで返ってきます。これは結局dict
みたいなものです。[*1]
軌道はy
に格納されていて、次元×STEP数のnp.ndarray
です。
3次元プロットをすることで軌道を確認します。
=
,
アトラクターを確認することができました。
method=
scipy.integrate.solve_ivp
ではOdeSolver
を使用して計算を行なっています。
solver
は引数method
で設定されます。デフォルトではRK45
が採用されています。
RK45はEmbedded Runge-Kutta Formulasの一つで、刻み幅の自動調整が行われています。
刻み幅h
の推移を見ることで調整しているかどうかを確認します。
刻み幅の動的更新が確認できます。
Embedded Runge-Kutta
Embedded Runge-Kutta(埋め込みルンゲ・クッタ)では $p$ 次の $b$ と、 $\hat{p}$ 次( $\hat{p}=p-1 \ or \ \hat{p}=p+1$ )の $\hat{b}$ を利用して計算を行います。
Butcher tableauは以下のようになります。
$$$ \def\arraystretch{1.5} \begin{array}{c|c} 0 & & & \\ c_1 & a_{21} & & & \\ c_2 & a_{31} & a_{32} & & \\ \vdots & \vdots & \vdots & \ddots & & \\ c_s & a_{s1} & a_{s2} & \cdots & a_{s, s-1} & \\ \hline & b_1 & b_2 & \cdots & b_{s-1} & b_s \\ \hline \\ & \hat{b_1} & \hat{b_2} & \cdots & \hat{b_{s-1}} & \hat{b_s} \\ \end{array} $$$
$ k_s =f(t_0+c_sh, y_0+h(a_{s1}k_1+...+a_{s, s-1}k_{s-1})) $
$p$次の$y_1$を計算します。
$ y_1 = y_0 + h(b_1k_1+...+b_sk_s) $
$\hat{p}$次の$\hat{y}_1$を計算します。
$ \hat{y}_1 = y_0 + h(\hat{b}_1k_1+...+\hat{b}_sk_s) $
solve_ivpはこの差 $y_1 - \hat{y}_1$ を$sc_i$以下に保ちます。
$ |y_{1i}-\hat{y_{1i}}|\leq sc_i \\ sc_i=Atol_i + max(|y_{0i}|,|y_{1i}|) Rtol_i $
rtol=1e-3, atol=1e-6
とデフォルトで設定されています。
# scipy/integrate/_ivp/rk.py
...
...
= + *
...
$Rtol$と$Atol$を使用して、相対許容誤差と絶対許容誤差を制御できる。らしいです。
確かに、$Rtol$と$Atol$を小さくすることで$sc_i$を小さくできるので、誤差を小さくすることは可能です。
ただ、値の設定に関してちゃんとした物を発見できていないのでそのうち追記します。
$sc_i$を使用した$err$を元に刻み幅h
を決定します。
$ err = \sqrt{\frac{1}{n}\displaystyle\sum_{i=1}^n({\frac{y_{1i}-\hat{y}_{1i}}{sc_i}})^2} $
RK45
C,A,B
に上記の記号$\bold{c}, \bold{a}, \bold{b}$が割り当てられています。E
は$\bold{\hat{b}}$に当たります。
こちらの値は以下の論文の値を参考にしていました。
Some Practical Runge-Kutta Formulas
...
=
=
=
=
...
刻み幅を指定する
刻み幅を指定する場合はt_eval
を渡します。
例えば刻み幅0.01で計算を行いたい場合は以下のようになります。
(np.linspace
を使用していますが、np.arange(0.0, 40.01, 0.01)
でも大丈夫です。)
=
実装を簡単に確認します。
(1) まずt_eval
を指定することでnp.searchsorted
を使用してt_eval
からt_eval_step
として抽出します。
...
= None
= # <- (1)
=
...
=
# <- (2)
=
(2) t_eval_step
とP
を使用して補間sol(t_eval_step)
*2が行われます。
これによりt_evalにおける軌道が計算されます。
(補間についてはそのうち調べます。)
参考文献
Solving Ordinary Differential Equations I
日本語訳も存在しますが、すでに販売を終了しており、図書館の貸し出し状態は不明だったのでこれを参考にしています。
II.1 THe First Runge-Kutta Methods
II.4 Practical Error Estimation and Step Size Selection
scipy.integrate.solve_ivp
scipy.integrate.solve_ivp | SciPy documentation
おまけ
OdeResult
scipy.solve_ivp
の返り値であるOdeResult
を確認します。
# scipy/integrate/_ivp/ivp.py
pass
# scipy/optimize/optimize.py
""" Represents the optimization result.
...
"""
return
=
=
= + 1
return
return +
return
sol(t_eval_step)
RkDenseOutput(...)(t_eval_step)
は結局RkDenseOutput(...)._call_impl(t_eval_step)
です。
# scipy/integrate/_ivp/rk.py
...
= # <-
...
=
return
...
# <- 補間処理部分
= /
=
=
=
=
= *
+=
+=
return
# scipy/integrate/_ivp/base.py
...
...
return # <- 結局ここ
...