Rパッケージ deSolveを使う

ここではRパッケージ deSolveの簡単な使い方を紹介する。詳しくはPackage deSolve: Solving Initial Value Differential Equations in R を参照。

対象とする常微分方程式系

$(x,y,z)\in \mathbb{R}^3$内の次の微分方程式（Lorenz系）を例にしてDeSolveを使ってみよう。 \begin{align} \frac{dx}{dt} &=\sigma (y-x)\\ \frac{dy}{dt} &=r x -y - xz\\ \frac{dz}{dt} &=-b z + xy \end{align} [Deterministic Nonperiodic Flow (Edward N. Lorenz, J. Atmos. Sci., 20(11963), 130–141）]。系の状態変数（state variables）は $(x, y, z)\in \mathbb{R}^3$、微分方程式系を定めるパラメータは$\sigma, b, r$ であるが、Lorenzに倣って $\sigma = 10$, $b=8/3$ を固定し、$r$ を様々に変化させるパラメータとしよう。微分方程式の初期状態（初期値 initial values）を $x(0)=0, y(0)= 0.1, z(0)= 0$ とする。

deSolveの微分方程式の取り扱い

パッケージ deSolveを使うためには、まず、パッケージ deSolve を次のようにして読み込む。以下で、Rのプロンプトを表す記号(>)は省略している。

> library(deSolve)

次いで、数値計算する対象の微分方程式を、パラメータと初期状態の組および解くべき時間幅とその差分刻みを定め、そして方程式を定義する。

パラメータと初期状態の組

Lorentz系でのパラメータを $\sigma(s) = 10$, $b=8/3 , r=28$ とするとき、これらをベクトル的に並べて任意のパラメータ変数名（ここでは parameters とした）を付ける。

parameters <- c(s = 10, b = 8/3, r = 28)

Lorenz系の初期条件値をベクトル的に並べて任意の変数名（ここでは initial とした）を付ける。

initial <- c(x=0, y=0.1, z=0)

解くべき時間幅とその差分刻み

微分方程式を解くべき時間（通常は時刻 $t=0$ から適当な時刻まで）をさだめ変数名（ここでは t とした）を付けた。ここでは、時刻 0 から 100 までとし、数値計算のための時間刻み $\Delta t=0.01$ とした。

times <- seq(0, 100, 0.01)

微分方程式の定義

deSolveでは微分方程式を次のようにして関数定義し、任意の名称（こｋでは Lorenz）を付けた。

Lorenz <- function(t, state, parameters) {
	with(as.list(c(state, parameters)), {
		dx <- s * (y - x)
		dy <- r * x - y - x * z
		dz <- -b * z + x * y
		list(c(dx, dy, dz))
	})
}

ここでの引数 t, state, parameters は、それぞれ「時間幅とその差分刻み」、「初期状態の組」および「パラメータの組」であるが、先に名付けた名称と同じである必要はないことに注意。先に名付けた変数名は、deSolveの数値積分ルーチンの呼び出し時に引数として与える（後述）。

deSolveで数値積分を実行する

以上の準備をして、次のようにして数値積分を実行する。その出力を out に代入している。

out <- ode(y = initial, times = times, func = Lorenz, parms = parameters)

ここで、引数値として initial, times, Lorenz, parameters は先に定めた変数値（任意の名称でよかった）を代入するのであるが、それぞれ、その代入先である y, times, func, parms は『このまま』の名称として使わねばならない。

以上まとめると、次のようなソースとなる。

library(deSolve)# 微分方程式の数値解法パッケージ

parameters <- c(s = 10, b = 8/3, r = 28)# パラメータのセット
initial <- c(x=0, y=0.1, z=0)# 初期条件
times <- seq(0, 100, 0.01)# 差分刻み

Lorenz <- function(t, state, parameters) {
	with(as.list(c(state, parameters)), {
		dx <- s * (y - x)
		dy <- r * x - y - x * z
		dz <- -b * z + x * y
		list(c(dx, dy, dz))
	})
}

# deSolveを使った数値積分を out に出力
out <- ode(y = initial, times = times, func = Lorenz, parms = parameters)

計算結果

出力された数値積分結果 out の構造を見てみよう

> str(out)
 deSolve [1:10001, 1:4] 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr [1:4] "time" "x" "y" "z"
 - attr(*, "istate")= int [1:21] 2 11134 22324 NA 7 7 0 68 23 NA ...
 - attr(*, "rstate")= num [1:5] 0.01 0.01 100.01 0 72.92
 - attr(*, "lengthvar")= int 3
 - attr(*, "class")= chr [1:2] "deSolve" "matrix"
 - attr(*, "type")= chr "lsoda"

計算結果はデータフレームとして書き出されており、先頭の6業は次のようである。

> head(out)
     time           x         y            z
[1,] 0.00 0.000000000 0.1000000 0.000000e+00
[2,] 0.01 0.009512481 0.1003535 4.791048e-06
[3,] 0.02 0.018277226 0.1032424 1.869053e-05
[4,] 0.03 0.026597741 0.1084756 4.168197e-05
[5,] 0.04 0.034733515 0.1159428 7.456224e-05
[6,] 0.05 0.042911665 0.1256036 1.188774e-04

一列目が経過時間、2列目以降が対応する時刻の状態変数値であることがわかる。