Python NumPy SciPy : 1 階常微分方程式の解法

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

関数の種類

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

Table 1: 常微分方程式の積分器
関数 説明
odeint ODEPACK の LSODE に基づいて計算する、シンプルなインターフェースの積分器
ode より一般的なインターフェースの積分器
complex_ode 複素方程式系用の ode のラッパー
Table 2: ode の数値解法
関数 数値解法 問題のタイプ
vode(method="adams") Adams 法(予測子修正子法) Nonstiff
vode(method="bdf") BDF (線型多段法) Stiff
zvode(method="adams") 複素方程式系用、Adams 法 Nonstiff
zvode(method="bdf") 複素方程式系用、BDF Stiff
lsoda Adams 法と BDF を動的に切り変える Stiff, Nonstiff
dopri5 5(4)次の Runge-Kutta 法(Dorman-Prince 法) Nonstiff
dop853 8(5,3)次の Runge-Kutta 法(Dorman-Prince 法) Nonstiff

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 次元プロットしてみます。
以下の数式を

150617lorenz.png

関数オブジェクトとして定義し、odeint に初期値やパラメータを引数に設定するだけです。

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


def func(v, t, p, r, b):
    return [-p*v[0]+p*v[1], -v[0]*v[2]+r*v[0]-v[1], v[0]*v[1]-b*v[2]]

p = 10
r = 28
b = 8/3
v0 = [0.1, 0.1, 0.1]
t = np.arange(0, 100, 0.01)

v = odeint(func, v0, t, args=(p, r, b))

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(v[:, 0], v[:, 1], v[:, 2])
# plt.show()

15061701.png

コメント

Comments powered by Disqus
書籍更新情報
2016-10-21
Pythonによる科学技術計算 基礎編
PDF版の販売を開始しました。
販売ページはこちら

2016-09-09
Pythonによる科学技術計算 基礎編
1.2版への更新が可能になりました。
サポートページはこちら
電子書籍
Pythonによる科学技術計算 基礎編
Kindle ストア、Leanpubで販売中です
Pythonによる科学技術計算 基礎編
PDF版の販売はこちら
同人誌
技術書典(2016/6/25)
Emacs/org-modeのPDF作成術
電子版をBOOTHで販売中です
Emacs/org-modeのPDF作成術
Share