Webページでは数式を表すためにLaTeX表記が使えるMathJaxを利用しています。 WebブラウザにはSafari/Chrome/Firefoxを使って下さい（IEでは表示できないようです。）

SciPyで微分方程式を解く

モジュール SciPyの数値積分には微分方程式の数値解法 scipy.integrate.odeint が含まれている。 FORTRANライブラリのodepackにあるlsodaを使っているので高速な計算が可能である。『硬い方程式』(stiff equation)にも対応しており汎用性が高い。詳しい使い方はSciPyのドキュメントscipy.integrate.odeintを参照。

問題とする微分方程式系は$n$変数の1階微分方程式と$n$個の初期値の組とする。 scipy.integrate.odeintの基本的な使い方は次のようである。

import scipy.integrate.odeint

orbit = scipy.integrate.odeint(vectorfield, y0, t, args)

引数 vectorfield は微分方程式の右辺を定めるベクトル場リスト（以下の例）、u0 は初期値のリスト（array）、tは解くべき時間点のarrayリスト（最初の時間点での変数値が初期値 u0）、args はvectorfield に渡す引数のタプルである。返り値 orbit はarrayリストで、解かれる時間点に応じた軌道変数値のarrayリスト [u0, u1, u2,...] である（各要素もまた変数の数の長さからなるarrayである）。

具体例

ベクトル場関数の与え方

Lorenz系（Deterministic Nonperiodic Flow J. Atmos. Sci., 20(11963), 130–141)では、odeintで呼ばれるベクトル場は次のように書くことができる。 p, b, r はベクトル場のパラメータ（ここでは定数）である。

    def lorenz(u, t, p, b, r):
        x = u[0]
        y = u[1]
        z = u[2]
        dxdt = -p * x + p * y
        dydt = -x * z + r * x - y
        dzdt = x * y - b * z
        return([dxdt, dydt, dzdt])

初期値の設定

odeintに渡す初期値は、$x_0 = 0.1, y_0=0.1, z_0=0.1$とすると次のように与える。

u0 = [0.1, 0.1, 0.1]

解くべき時間点

初期値 u0 の時間点を 0.0 としたとき、差分刻み 0.01 で時刻 40 までの計算のために odeint に渡す時間点 times は次のように与える。

dt = 0.01
T = 40.0
times = np.arange(0.0, T, dt)

プログラム例

以上をまとめて、Lorenz系をSciPyを使って微分方程式の数値解を求めて、軌道要素を3次元プロットするプログラム scipy_diff.py を次のように書く。 odeintで返される軌道要素は多数行$\times 3$のarrayであり、Numpyを使って効率化を図るで取り上げたように、各列全ての要素リストはそれぞれ orbit[:, 0], orbit[:, 1], orbit[:, 2] でスライスできることに注意しよう（コメントアウトされた xorbit, yorbit, zorbitを改めて抜き出す計算は必要はない）。

# -*- coding: utf-8 -*-
# numerila solution using scipy.integrate
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D

def lorenz(u, t, p, b, r):
    x = u[0]
    y = u[1]
    z = u[2]
    dxdt = -p * x + p * y
    dydt = -x * z + r * x - y
    dzdt = x * y - b * z
    return([dxdt, dydt, dzdt])

# parameters
p =10
b = 8.0 / 3.0
r = 28
# initial values
u0 = [0.1, 0.1, 0.1]
# setting
dt = 0.01
T = 40.0
times = np.arange(0.0, T, dt)
# numerical solution using scipy.integrate
args = (p, b, r)
orbit = odeint(lorenz, x0, times, args)

#xorbit = []
#yorbit = []
#zorbit = []
#for x, y, z in orbit:
#    xorbit.append(x)
#    yorbit.append(y)
#    zorbit.append(z)

#plot3D using matplotlib
fig = plt.figure()
ax = fig.gca(projection='3d')
#ax.plot(xorbit, yorbit, zorbit)
ax.plot(orbit[:, 0], orbit[:, 1], orbit[:, 2])
ax.set_xlabel('x(t)')
ax.set_ylabel('y(t)')
ax.set_zlabel('z(t)')
plt.show()