6/13　MIKUセミナー

R 数理モデル 数学

MikuHatsune2012-06-13

数学セミナーの10月号「突発的変化の数理」。
　
以前、人口問題を扱ったが、ここでも生物の個体数の時間変化を記述したロジスティク方程式
f(N)=r(1-\frac{N}{k})
を扱っている。
1978年に、害虫の突発的大発生モデルとして
N'(t)=r(1-\frac{N}{k}) - \frac{\beta N^2}{\alpha ^ 2+N^2}
\frac{\beta N^2}{\alpha ^ 2+N^2}：害虫を餌とする捕食者による害虫駆除の効果
r>0：害虫の増殖率（出生率と死亡率の差）
k>0：可食され得る葉っぱの量に関連した扶養能力
\alpha：飽和状態が始まるときの害虫の密度の大きさ
\beta：ゆっくりとした森の変化に依存した値
というものがあるらしい。簡単に
N'(t)=r(1-\frac{N}{k}) - \frac{N^2}{1+N^2}
として、N'(t)=0となる平衡点を考える。
r(1-\frac{N}{k}) = \frac{\beta N^2}{\alpha ^ 2+N^2}
をガリガリNについて解くと、
g(N)=rN^3-rkN^2+(k+r)N-rk
r、kを色々動かした時に、g(N)=0となるNがいくつあるかを愚直に数えてみる。

#N, k, rをそれぞれ動かして計算する。
a1 <- function(N, k, r){
	r*N^3 - r*k*N^2 + (k+r)*N -r*k 
}
N <- seq(0, 100, length=3000) #理想としては∞
k <- seq(0, 10, length=300) #大きければ大きいほどいいかもしれない
r <- seq(0, 1, length=100) #実験的に0~1まで与えればいいことがわかったので
A <- matrix(0, length(k), length(r)) #データ格納

for(k1 in 1:length(k)){
for(r1 in 1:length(r)){
	w <- a1(N, k[k1], r[r1])
	#g(N)=0の判定には、中間値の定理を用いる。
	#g(i), g(i+1)の積が負なら、必ず0を越さないといけないというあれ。
	A[k1, r1] <- sum(head(w, -1L) * tail(w, -1L) < 0)
}
}
#Nが0個なら白、1個なら緑、2個なら青、3個なら赤。
cols <- c("white", "green", "blue", "red")
image(k, r, A, col=cols)
#3次方程式なので最大個数は3つのはず。
range(A)


基本的にはNの個数は1個か3個。
イメージ的には、ふにゃふにゃした紙を上から見下ろして、視線が平面を何回通過するかを数えている感じ。