今回も書籍『A Biologist's Guide to Mathematical Modeling in Ecology and Evolution』で学んだ内容を復習します。内容としては、計算モデルの収束判定です。

捕食者-非捕食者モデル

今回例とするのは捕食者-非捕食者モデルになります。環境中の捕食者（キツネなど）と非捕食者（ウサギなど）の個体数の増減を示すモデルになります。

\[ \frac{dn_1}{dt}=\theta - acn_{1}(t)n_{2}(t) \\ \frac{dn_2}{dt}=\epsilon acn_{1}(t)n_{2}(t) - \delta n_{2}(t) \\ \]

また、各記号の意味は↓になります。

記号	意味
$n_{1}(t)$	非捕食者の個体数
$n_{2}(t)$	捕食者の個体数
$\theta$	非捕食者の流入速度
$a$	捕食が発生する確率
$c$	捕食者と非捕食者が接触する速度
$\epsilon$	捕食数を捕食者の個体数増に変換する係数
$\delta$	捕食者の死亡速度

ベクトル場の記述

まず、ベクトル場を書いてみます。各パラメーターの値は以下にしています。

記号	意味	値
$\theta$	非捕食者の流入速度	10
$a$	捕食が発生する確率	1
$c$	捕食者と非捕食者が接触する速度	0.1
$\epsilon$	捕食数を捕食者の個体数増に変換する係数	1
$\delta$	捕食者の死亡速度	0.3

ENV["GKS_ENCODING"] = "utf-8"
using Plots
gr()
δ=0.3;
θ=10;
ϵ=1;
ac=0.1;

nextn1(n1,n2,θ,ac)= θ - ac*n1*n2;
nextn2(n1,n2,δ,ϵ,ac)= ϵ*ac*n1*n2 - δ*n2;

N = zeros(10,10);
global p1 =plot(xlabel="n1",ylabel="n2",xlims=(0,45),ylims=(0,45));

for i in 1:20 ,j in 1:20
  n1 = 2*i;
  n2 = 2*j;
  scale = 100;
global p1=plot!([n1,n1+nextn1(n1,n2,θ,ac)/scale],[n2,n2+nextn2(n1,n2,δ,ϵ,ac)/scale],
      arrow = arrow(0.4,0.4), label="",linecolor = "black")
end
p1

この条件では、安定的な平衡点が存在しそうですね。

微分方程式の解

次に、この微分方程式を解いてみます。

using DifferentialEquations

function predatorprey(dn,n,p,t)
  dn[1]= p[1] - p[2]*n[1]*n[2];
  dn[2]= p[3]*p[2]*n[1]*n[2] - p[4]*n[2];
end

n0 = [1.0;1.0];
tspan = (0.0,40.0);
p=[θ,ac,ϵ,δ];
prob = ODEProblem(predatorprey,n0,tspan,p)
sol =solve(prob)
lastpoint = sol(40)
p1 = plot!(sol,vars=(1,2),xlims=(0,45),ylims=(0,45),linewidth=4,linecolor = :red,label="")
p1 = scatter!([lastpoint[1]],[lastpoint[2]],color = :red, markersize = 7,label = "")

p0 = plot(sol,vars=(0,1),label = "n1")
p0 = plot!(sol,vars=(0,2),label = "n2",xlabel="Time", ylabel="Population size")

p0

時系列データを見ると値が収束していっていることがわかります。

収束の様子をベクトル場に重ねてみますと、矢印の流れに乗って収束していっていることがわかります。

一方、ここで$δ=2$に変更するとグラフの様子が様変わりします。

パラメーターを$δ=2$に変更したことにより。時系列データの値が振動しながら収束していっています。

また、ベクトル場にプロットすると渦を巻いているようになります。

収束の判定方法

パラメーターの値の違いによる収束への影響を調べてみます。今回の式を再掲します。

\[ f_1 = \frac{dn_1}{dt}=\theta - acn_{1}(t)n_{2}(t) \\ f_2 = \frac{dn_2}{dt}=\epsilon acn_{1}(t)n_{2}(t) - \delta n_{2}(t) \\ \]

平衡点を求める

　平衡点への収束の様子を調べるには、まず平衡点$(\hat{n_1},\hat{n_2})$を求める必要があります。ここで平衡点とは$\frac{dn_1}{dt} = 0$ かつ $\frac{dn_2}{dt} = 0$となる点です。

\[ 0=\theta - ac\hat{n_1}\hat{n_2} \\ 0=\epsilon ac\hat{n_1}\hat{n_2} - \delta \hat{n_2} \\ \]

このとき、$\hat{n_2} \neq 0$とすると平衡点は以下となります。

\[ \hat{n_1} = \frac{\delta}{\epsilon ac}\\ \hat{n_2} = \frac{\theta \epsilon}{\delta} \]

ヤコビ行列を求める

つづいて、この平衡点$(\hat{n_1},\hat{n_2})$ におけるヤコビ行列を求めます。ヤコビ行列とは以下のものを指します。

\[ J = \left( \begin{array}{cc} \frac{\partial f_1}{\partial n_1} & \frac{\partial f_1}{\partial n_2} \\ \frac{\partial f_2}{\partial n_1} & \frac{\partial f_2}{\partial n_2} \\ \end{array} \right)\\ \]

したがって今回の場合のヤコビ行列は以下になります。

\[ J = \left( \begin{array}{cc} -ac\hat{n_2} & -ac\hat{n_1} \\ \epsilon ac\hat{n_2} & \epsilon ac\hat{n_1} - \delta \\ \end{array} \right)\\ \]

ここで$\hat{n_1} = \frac{\delta}{\epsilon ac}, \hat{n_2} = \frac{\theta \epsilon}{\delta}$のため上のヤコビ行列は以下のように変形できます。

\[ J = \left( \begin{array}{cc} -ac\epsilon \theta / \delta & -\delta / \epsilon \\ \epsilon^2 ac\theta /\delta & 0 \\ \end{array} \right)\\ \]

収束判定

長くなりましたが、このヤコビ行列を用いて収束の様子を判定することが出来ます。

ヤコビ行列の全ての固有値の実部が負のとき、安定収束
上記の場合でも、固有値に 複素数 が含まれている場合、時系列で振動し、ベクトル場では渦を巻く

では$\delta = 0.3$および$\delta = 2$の場合のヤコビ行列の固有値を求めてみます。

まず$\delta = 0.3$の場合を見てみます。

using LinearAlgebra
δ=0.3;
θ=10;
ϵ=1;
ac=0.1;

J = [-ac*ϵ*θ/δ -δ/ϵ; ϵ^2*ac*θ/δ 0 ];
eigvals(J)

2-element Array{Float64,1}:
 -3.0               
 -0.3333333333333333

計算の結果$\delta = 0.3$の時の固有値は全て負の実数なので、確かに安定的に収束することがわかります。

また、$\delta = 2$の場合の固有値は以下のようになります。

using LinearAlgebra
δ=2;
θ=10;
ϵ=1;
ac=0.1;

J = [-ac*ϵ*θ/δ -δ/ϵ; ϵ^2*ac*θ/δ 0 ];
eigvals(J)

2-element Array{Complex{Float64},1}:
 -0.24999999999999997 - 0.9682458365518543im
 -0.24999999999999997 + 0.9682458365518543im

上記のように$\delta = 2$の場合は固有値は全て実部が負の複素数です。なので、安定的に収束するものの振動や渦巻が確認できます。

では、上記の場合はどうでしょうか。$\epsilon = -1$としました。

using LinearAlgebra
δ=2;
θ=10;
ϵ=-1;
ac=0.1;

J = [-ac*ϵ*θ/δ -δ/ϵ; ϵ^2*ac*θ/δ 0 ];
eigvals(J)

2-element Array{Float64,1}:
 -0.7807764064044151
  1.2807764064044151

この場合、実部が正の固有値が含まれることになります。

この場合、収束せず発散することになります。

今回はここまでにします。

システムとモデリング

modelica, Julia, Design Structure Matrix, SysML, 他モデリング全般について。

【Julia入門】Juliaで計算モデルの解析：収束判定