概要

JuliaにはSundials.jlという微分方程式の数値計算のためのライブラリがあります。

qiita.com

Sundialsに関しては上のリンクにかかれていますが、

SundialsライブラリはJulia専用ではなく、いろんな言語で使用できて、Matlabでも使うことができる

ような微分方程式のためのライブラリのようです。

今回はこのSundials.jlを使ってロトカ・ヴォルテラ方程式の数値計算をしてみたいと思います。

実験

これからJuliaで数値計算の実験をしていきたいと思いますが、まずはロトカ・ヴォルテラ方程式に関して簡単に見ていきたいと思います。

ロトカ・ヴォルテラ方程式とは

今回のテーマとして設定したロトカ・ヴォルテラ方程式とは何かと言うとWikipediaによると、

生物の捕食-被食関係による個体数の変動を表現する数理モデルの一種

というものです。

ロトカ・ヴォルテラの方程式 - Wikipedia

食うものと食われるものがいた時に、それぞれの個体数が時間によってどのように変化していくかをシミュレーションする時に使えそうな数理モデルということでしょうか*1。

このモデルに基づいて可視化をした一つの例がWikipediaに乗っていますが、

f:id:rf00:20181124112928p:plain

この図のように、横軸に時間、縦軸に個体数をとって、青が食うもの赤が食われるものとした時に、食うものが増えると食われるものが減り、その状態から段々と食うものが減り食われるものが増え食うものが減る、そしてそれが交互に繰り返されていく、といった現象をモデル化している事がわかります。

このモデルによって表現できそうなものを考えると魚と漁師の関係なんか表現できるかもしれないですね。

そしてロトカ・ヴォルテラ方程式は以下のように定義されているものです。

${\frac{dx}{dt} = ax - bxy}\\ {\frac{dy}{dt} = cxy - dy}$

ここで x は被食者の個体数、 y は捕食者の個体数、t は時間をあらわし、4つの係数 a, b, c, d は正の実数のパラメータである。

というのがWikipediaからの説明です。

Juliaでの実践

それではJuliaでの実践に移ります。Sundialsライブラリは入っている前提で進めます。

まずはSundialsで計算するためにロトカ・ヴォルテラ方程式を関数として定義します。

using Sundials
using Plots 

#ロトカヴォルテラ方程式
function lotka_volterra(t, u, du)
    du[1] = u[1] * (a   - b * u[2]) # dx/dt
    du[2] = u[2] * (c * u[1] - d) # dy/dt
    du
end

若干式変形をしていますが定義どおりの実装です。

ハイパーパラメータの設定は以下の通りにします。

# パラメータ設定
a = 1.5
b = 1.0
c = 1.0
d = 3.0
u0 = [10.0; 5.0] #初期値
t = [0.0:0.1:10.0;] #時間

そしてSundialsで数値計算するには以下を実行します。

#数値計算
res = Sundials.cvode(lotka_volterra, u0, t)
res

そうすると以下のような結果が返ってきます。

101×2 Array{Float64,2}:
 10.0       5.0     
  5.94677   8.23517 
  2.79878   9.31826 
  1.32512   8.39552 
  0.720449  6.85492 
  0.454605  5.37652 
  0.327537  4.14037 
  0.264391  3.15862 
  0.232879  2.39858 
  0.219388  1.81738 
  0.217489  1.37597 
  0.22413   1.04198 
  0.237742  0.789809
  ⋮                 
  1.28688   8.28976 
  0.704551  6.75654 
  0.448128  5.29459 
  0.32622   4.07489 
  0.264547  3.10902 
  0.234443  2.36062 
  0.221485  1.78908 
  0.220224  1.35475 
  0.227247  1.02634 
  0.24143   0.778273
  0.262214  0.591056
  0.289233  0.450104

これをプロットすると、

# プロット
x = res[:,1]
y = res[:,2]

plot(t, x, label='x')
plot!(t, y, label='y')

png("lotka_volterra")

f:id:rf00:20181124115606p:plain

このようにいい感じに食われるものが最初は多いが食われていって減少、食うものが増加しするが、食われるものが減っていくため減少し、食われるものが増加するという関係性がモデル化できているようですね！

まとめ

今回はロトカ・ヴォルテラ方程式の数値計算をJuliaのSundialsライブラリを使ってやってみました。

とても簡単にできることが確認できたので個人的には今後微分方程式の理解を深めるために積極的に使ってみたいと思います。

それでは間違いなどありましたら指摘お願いいたします。

*1:個人的な思い出ではありますが大学院生の時に進化ゲームの文脈で聞いたことのあるものでもあります

Deep Karmaning

技術系の話から日常のことまで色々と書きます

JuliaのSundialsライブラリを使ってロトカ・ヴォルテラ方程式の数値計算をしてみた

概要

実験

ロトカ・ヴォルテラ方程式とは

Juliaでの実践

まとめ