Python Numpy Scipy : 1 階常微分方程式の解法

Scipy.integrate パッケージ の 1 階常微分方程式の積分器についてメモ。

多変量

スカラー関数 結果がスカラー

The specific optimization method interfaces below in this subsection are not recommended for use in new scripts; all of these methods are accessible via a newer, more consistent interface provided by the functions above. このサブセクションの下の特定の最適化方法インタフェースは、新しいスクリ プトで使用することはお勧めしません。これらの方法のすべては、上記の機能 が提供する新しい、より一貫性のあるインタフェースを介してアクセス可能で す。

最小二乗最小化(leastsq)とカーブフィッティング(curve_fit) スカラ単変量関数 minimizers(minimize_scalar)とルートファインダー(newton)

関数の種類

1 階常微分方程式を解く関数は表 1 のものがあります。基本的には odeint を使い、ode は set_integrator メソッドで様々なオプションを指定して表 2 の数値解法で計算したいときに使います。

Table 1: 多変量スカラ関数の制約なしと制約付き最小化(minimize)
関数 説明 制約条件 導関数
minimize(method="Anneal") Aneal 法 非制約  
minimize(method="dogleg") 信頼領域 dog-leg 法 非制約  
minimize(method="trust-ncg") Newton conjugate gradient trust-region algorithm 非制約  
minimize(method="L-BFGS-B") Limited-memory BFGS 制約  
minimize(method="TNC") truncated Newton algorithm 制約  
minimize(method="COBYLA") 線形近似による制約付き最適化 制約 不要
minimize(method="SLSQP") Sequential Least SQuares Programming 制約  
minimize(method="custom") custom    

Anneal 精度はよくないが、とける問題がある。

dogleg

trust-ncg

L-BFGS-B

TNC

COBYLA 線形計画みたい。FORTRAN。不等式の制約条件がある場合はこれを使う。非線形制約。

SLSQP

custom

Table 2: 多変量方程式系ソルバー(root)
関数 説明
fsolve MINPACK の HYBRD,HYBRJ。修正 powell
root(method="hybr") MINPACK の HYBRD,HYBRJ。修正 powell
root(method="lm") 修正 Levenberg-Marquard
root(method="broyden1") 非厳密ニュートン法 broyden 1st jacobian
root(method="broyden2") 非厳密ニュートン法 royden 2nd jacobian
root(method="anderson") 非厳密ニュートン法 nderson mixing
root(method="linearmixing") 非厳密ニュートン法 ninear mixing
root(method="diagbroyden") 非厳密ニュートン法 iagonal broyden
root(method="excitingmixing") 非厳密ニュートン法 exciting mixing
root(method="krylob") 非厳密ニュートン法 krylob

2 通り以上の方法を試して比較すること。

Table 3: Equation (Local) Minimizers
関数 説明
leastsq BDF (線型多段法)
nnls BDF (線型多段法)
Table 4: グローバル最適化ルーチン(aneeal, basinhopping, differential_evolution)
関数 説明
basinhoppping BDF (線型多段法)
brute BDF (線型多段法)
differential_evolution BDF (線型多段法)

aneal is deprecated

Table 5: グローバル最適化ルーチン(aneeal, basinhopping, differential_evolution)
関数 説明
rosen BDF (線型多段法)
rosen_der BDF (線型多段法)
rosen_hess BDF (線型多段法)
rosen_hess_prod BDF (線型多段法)
Table 6: curve_fitting
関数 説明
curve_fit BDF (線型多段法)
Table 7: 解探索スカラー関数
関数 説明
brentq BDF (線型多段法)
brenth BDF (線型多段法)
ridder BDF (線型多段法)
bisect BDF (線型多段法)
newton  
fixed_point  
Table 8: linear programming
関数 説明
inprog BDF (線型多段法)
Table 9: utility functions
関数 説明
approx_fprime BDF (線型多段法)
bracket  
check_grad  
line_search  
show_options  

odeint は odepack の LSODE ルーチン を使用しています。LSODE は問題が Nonstiff なら Adams 法(Adams-Bashforth-Moulton 法)、Stiff なら BDF(後退微分公式)で計算します。問題が Stiff かどうか気にする必要がありませんし、計算効率や計算精度も比較的高いため、まず最初は odeint で計算してみましょう。

vode はオプションの method に指定することで Adams 法と BDF を選んで計算できます。odeint で選択されなかった方法についても計算してみて結果を比較したりするときに使います。

lsoda は LSODE に似た LSODA ルーチンで計算します。最初は Nonstiff の Adams 法で計算を始め、 計算中に必要に応じて BDF に切り変えて計算します。そのため odeint と同じでどちらの問題のタイプも解くことができます。問題によって LSODE と LSODA で計算効率が異なるので、計算効率に拘るなら lsoda でも計算してみて結果を比較してみるとよいと思います

dopri5 は Matlab の ode45 にあたる、5(4)次の Runge-Kutta 法で計算します。odeint がうまくいかずその理由がわからないときや、平易な問題で計算効率を気にしないときなどに使えばいいと思います。

dop853 は Nonstiff の問題で、多少計算時間が長くても高精度な解を求めたいときに使用します。

数値計算法の名著 NUMERICAL RECIPES in C では初期ステップの計算の仕方をどうするか、刻み幅も次数も変化するようなプログラムを書くのは大変であることから、予測子修正子法は薦められない方法としていました。しかし、今ではその辺の問題も解消され、予測子修正子法は効率、精度共に高い方法として広く使われています。中身はブラックボックス的でもありがたく使わせてもらいましょう。

使い方の例

定番のローレンツ方程式を 3 次元プロットしてみます。以下の数式を

# coding: utf-8
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import j1
from scipy.optimize import minimize_scalar

res1 = minimize_scalar(j1, method='brent')
print(res1)

res1b = minimize_scalar(j1, bracket=(-25,-10,1), method='brent')
print(res1b)

res2 = minimize_scalar(j1, bounds=(4, 7), method='bounded')
print(res2)

x = np.arange(-30, 10, 0.1)

plt.figure()
plt.plot(x, j1(x))
plt.plot(res1.x, res1.fun, "ob", label="b")
plt.plot(res1b.x, res1b.fun, "oy", label="y")
plt.plot(res2.x, res2.fun, "or", label="r")
plt.legend()
plt.show()

15062101.png

コメント

Comments powered by Disqus