juliaで微分方程式を解くチュートリアルを書いてみました。
数理生物学に興味のある読者を想定して、ロトカ-ヴォルテラ方程式を題材に選んでいます。
- 生命現象を数理的に把握したい
- なんとなくかっこいいからやってみたい
そんな人に向けて書かれた記事です。
理論的背景はともかく、まずは julia を高性能な電卓として使ってみましょう。
tldr
完成品は「可読性が高い関数定義のやり方」に置いてあります。
数学的な導入
微分方程式
たとえば、いま興味のある量が2つくらいあり、そのどちらも時間に依存して変化するとします。
単なる表記の問題ですが、興味のある量 $x,y$ をひとまとめにしてベクトル量 $\boldsymbol {u} = (x,y) $ と書きましょう。
この $\boldsymbol {u} $ の時間変化率 $d\boldsymbol {u}/dt$ が、適当な関数 $f$ によって表されるとします。
\begin{align}
\dfrac {d}{dt} \boldsymbol {u} &= f\left( \boldsymbol {u},\boldsymbol {p},t\right) \\
\end{align}
ここでベクトル $\boldsymbol {p} $ は種々のパラメータ、変数 $t$ は時刻です。
任意の時刻 $t$ において上の式の等号が成立するとき、 $\boldsymbol {u} $ はどういった形をしているのでしょうか？
このような疑問を考えるとき、与えられた等式は微分方程式とみなすことができます。
ロトカ-ヴォルテラ方程式
ロトカ-ヴォルテラの捕食―被捕食モデルは、次の連立微分方程式で表されます。
\begin{align}
\frac{dN}{dt} &= -a NP + b N ~~&(1)\\
\frac{dP}{dt} &= c NP - d P ~~&(2)\\
\end{align}
ここで、$N$、$P$ はそれぞれ被捕食者の個体数、捕食者の個体数を表します。
$N$ は食物 nutrition の略、$P$ は捕食者 predator の略といったところでしょうか。
積 $NP$ は被捕食者と捕食者の遭遇回数を意味していると考えることができます。
被捕食者 N は遭遇回数が多ければ食べられてしまう危険が増すので、 $NP$ に比例して個体数が減少します。
これを表したのが式$(1)$の第一項 $-aNP$ です。
他方で、捕食者 P は捕食によって栄養を得、繁殖の機会が増えるので、個体数は遭遇回数に比例して増加します。式$(2)$第一項はこれを表し、個体数の変化は比例定数を $c$ として $+cNP$ です。
式$(1)$第二項 $+bN$ は被捕食者の自然発生率を、式$(2)$第二項 $-dP$ は捕食者の自然死の率を意味します。
プログラムのための数式へと書き換える
ロトカ-ヴォルテラ方程式を $ {u'} = f\left( {u}, {p},t\right)$ という形式に書き直してみましょう。
まず、興味のある量は $N$ と $P$ であり、この2つの変数をひとまとめにして $\boldsymbol {u} = (N,P) $ と書くのでした。
また、それ以外に式$(1)$、式$(2)$に含まれている $\boldsymbol {p} = (a, b, c, d)$ がパラメータとなります。
つまり、
\begin{align}
\boldsymbol {du} &:=\dfrac {d}{dt} \boldsymbol {u}= \left( \dfrac {dN}{dt},\dfrac {dP}{dt}\right) \\
\boldsymbol {u} &\,= \left( N,P\right) \\
\boldsymbol {p} &\,= \left( a,b,c,d\right) \\
\end{align}
となります。
さらに、これらの定義と式$(1)$と式$(2)$とを見比べて、スカラー量 $u,;p,;t$ の式で表していきます。
いま便宜的に、添字に $i$ と書くことを、ベクトルの $i$ 番目の要素の取り出し操作とみなします。
つまり、 ${du_{1}}=dN/dt,; {du_{2}}=dP/dt,; {u_{1}}=N,; {u_{2}}=P,; p_{1}=a,; p_{2}=b,; p_{3}=c,; p_{4}=d;$ です。
このことに注意して式$(1)$と式$(2)$を書き直したのが、以下の式$(1)'(2)'$です。
\begin{align}
du_{1} &= -p_{1}u_{1}u_{2}+p_{2}u _{1} ~~&(1)'\\
du_{2} &= p_{3}u_{1}u_{2}-p_{4}u _{2} ~~&(2)'\\
\end{align}
見比べてみると、全く同じ形になっていることが確認できますね。
実際に方程式を解いてみよう
先のように書き換えた微分方程式を julia の関数として定義します。
julia> # Lotka-Volterra 方程式の定義
julia> function lotka_volterra!(du,u,p,t)
du[1] = - p[1]*u[1]*u[2] + p[2]*u[1]
du[2] = p[3]*u[1]*u[2] - p[4]*u[2]
end
lotka_volterra! (generic function with 1 method)
初期条件とパラメータの指定
初期条件として、時刻 $t=0$ のときの被捕食者、捕食者それぞれの個体数 $\boldsymbol {u}(0) = \boldsymbol u_{0}$ を与えます。
個体数は整数ですが、小数点をつけて明示的に Float型で定義してやるのが julia においては大事です。
\begin{align}
\boldsymbol u_{0} = \boldsymbol {u}(0) = \left( N_{0}, P_{0}\right)= \left( 10, 7\right) \\
\end{align}
また、今回はパラメータ $\boldsymbol {p} = (a,b,c,d)$ を以下のように指定しました。
\begin{align}
\boldsymbol {p} = \left( a,b,c,d \right)= \left( 0.3,2.0,0.4,3.0\right) \\
\end{align}
以下では初期値 $\boldsymbol u_{0},;$ パラメータ $\boldsymbol {p}$ とシミュレート時間 $t = [0, 10]$ を定義しています。
julia> # 初期値 N0 = 10.0, P0 = 7.0 を与える
julia> # Float64型にするための 「.0」 が大切
julia> u0 = [10.0, 7.0]
2-element Array{Float64,1}:
10.0
7.0
julia> # パラメータ a, b, c, d を与える
julia> p = [0.3, 2.0, 0.4, 3.0]
4-element Array{Float64,1}:
0.3
2.0
0.4
3.0
julia> # シミュレートする区間 tspan の定義
julia> tspan = (0.0, 10.0)
(0.0, 10.0)
パッケージのインストールと書式
REPLで]と入力するとパッケージ管理モードに入ります。
微分方程式を解くためのライブラリ DifferentialEquations を追加します。
(v1.0) pkg> add DifferentialEquations
(v1.0) pkg> add Plots
(v1.0) pkg> precompile
DifferentialEquationsを使って微分方程式を解く際、書くのはたったの二行だけ（！）です。
julia> prob = ODEProblem(f,u0,tspan,p)
julia> sol = solve(prob)
解オブジェクト sol は解についての完全な情報を持っているそうです。
詳しくは公式ドキュメントを参照してください。
微分方程式を解いてグラフ描画まで
やっていきましょう。
なお、plot() 使用時に $ sudo apt install qt5-default（Debianの場合）が必要かもしれません。
julia> using DifferentialEquations
julia> using Plots;gr()
julia> problem = ODEProblem(lotka_volterra!,u0,tspan,p)
ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true
timespan: (0.0, 10.0)
u0: [10.0, 7.0]
julia> sol = solve(problem)
retcode: Success
Interpolation: Automatic order switching interpolation
t: 29-element Array{Float64,1}:
...（省略）...
u: 29-element Array{Array{Float64,1},1}:
...（省略）...
julia> plot(sol)
このグラフは縦軸が個体数、横軸は時間です。
個体数 $N,; P$ いずれも同じ周期で振動しているのがわかりますね。
また、 $u_{2} = P$ の位相が $u_{1} = N$ より遅れていることも確認できます。
解が $N-P,$ 平面の上でどういった図形を描くかを調べることもできます。
julia> plot(sol,vars=(1,2))
石ころみたいですね。
時間変化によって上の図形が描かれる過程をアニメーションで見ることもできます。
以下ではImageMagick を追加でインストールし、gifアニメを描画しています。
(v1.0) pkg> add ImageMagick
julia> using ImageMagick
julia> animate(sol,vars=(1,2),fps=30,lw=2,every=1)
今日のところはここまでとします。
可読性が高い関数定義のやり方
本稿の「プログラムのための数式へと書き換える」でやったことは実は少し冗長でした。
An Intro to DifferentialEquations.jl には以下のような文が出てきます。
However, that can be hard to follow since there's a lot of "programming" getting in the way.
先ほど「プログラムのための書き換え」でやったような小手先の操作は、できるなら避けたいですよね。
追加ライブラリの ParameterizedFunctions は、この問題を解決することができます。
このライブラリにおいて定義された@ode_defマクロを利用することで、関数を遥かにわかりやすく定義することができ、可読性の向上につながります。
下のコードがマクロ@ode_defを利用した関数定義の例です。
(v1.0) pkg> add ParameterizedFunctions
(v1.0) pkg> precompile
julia> using ParameterizedFunctions
julia> lv! = @ode_def LotkaVolterra begin
dN = - a*N*P + b*N
dP = c*N*P - d*P
end a b c d
(::LotkaVolterra) (generic function with 2 methods)
数式と見比べてもほぼ一緒ですね。これはすごい。
\begin{align}
\frac{dN}{dt} &= -a NP + b N ~~&(1)\\
\frac{dP}{dt} &= c NP - d P ~~&(2)\\
\end{align}
この定義でも solve() がバッチリ解いてくれます。
以下では先ほどと同様の手順を踏んで、同じ画像が出力されることを確認しています。
julia> using DifferentialEquations
julia> using ImageMagick
julia> u0 = [10.0, 7.0]
julia> p = [0.3, 2.0, 0.4, 3.0]
julia> tspan = (0.0, 10.0)
julia> problem = ODEProblem(lv!,u0,tspan,p)
julia> sol = solve(problem)
julia> animate(sol,vars=(1,2),fps=30,lw=2,every=1)
以上です。
あとがき
この記事の内容のほとんどは、公式ドキュメントなどの優れた英語版 Introduction から拾ってきたものです。
（特に [Intro to solving differential equations in Julia]
(https://www.youtube.com/watch?v=KPEqYtEd-zY&feature=youtu.behttps://qiita.com/tomotagwork/items/9e86d2fc67a7c7b644d5) と
An Intro to DifferentialEquations.jl は大いに参考にさせて頂きました）
下の参考サイトからは、本内容のすべてに加えて、もっと発展的な内容を学ぶことができます。
関連図書
『演習で学ぶ生命科学』第2版 羊土社（8章 マクロスケールのダイナミクス）
参考サイト
Intro to solving differential equations in Julia
An Intro to DifferentialEquations.jl
DifferentialEquations.jl 公式ドキュメント - Ordinary Differential Equations
DifferentialEquations.jl 公式ドキュメント - Solution Handling