ルードヴィッヒモデル

数理モデル

検索してもあまりいいページ出てこなかったですが式は以下の通り
\dot{n}=rn(1-\frac{n}{k})-\frac{sn^{2}}{1+n^{2}}

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

#n:個体数 k:飽和数 s:食傷捕食数 r:再生産率
parameters <- c(k = 100, s = 1, r=0.1)# パラメータのセット
initial <- c(n=5)# 初期条件
times <- seq(0, 100, 1)# 差分刻み

Ludwig <- function(t,state, parameters) {
	with(as.list(c(state, parameters)), {
		dn <- r*n*(1-n/k)-s*n^2/(1+n^2)
		return(list(c(dn)))
	})
}

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

parameters <- c(k = 100, s = 1, r=0.2)
out2 <- ode(y = initial, times = times, func = Ludwig, parms = parameters)
parameters <- c(k = 100, s = 1, r=0.3)
out3 <- ode(y = initial, times = times, func = Ludwig, parms = parameters)
parameters <- c(k = 100, s = 1, r=0.4)
out4 <- ode(y = initial, times = times, func = Ludwig, parms = parameters)
parameters <- c(k = 100, s = 1, r=0.5)
out5 <- ode(y = initial, times = times, func = Ludwig, parms = parameters)

out<-cbind(out1,out2[,2],out3[,2],out4[,2],out5[,2])

matplot(out,type="l")

r=0.202前後でカタストロフィー変化が生じる
このモデルでは初期状態n_{0}と再生産率rで\dot{n}の値が決まる