廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

deSolve パッケージでロトカ・ヴォルテラの方程式を数値的に解く

捕食者と被食者の増減関係

生物の種は互いに影響しあう.

ここでは食う‐食われる(捕食者と被食者;predator-prey)の関係を考える.

  • 被食者は捕食者がいなければ, その個体数に比例して増加する( dx/dt= \alpha x
  • 捕食者は被食者がいなければ, その個体数に比例して減少する( dy/dt= - \gamma y
  • 被食者は捕食されることによって, 自身の個体数および被食者の個体数に比例して減少する( -\beta x y
  • 捕食者は捕食によって, 自身の個体数および被食者の個体数に比例して増加する( +\delta xy

としてモデル化すると,

 \displaystyle \frac{dx}{dt} = \alpha x - \beta x y

 \displaystyle \frac{dy}{dt} = -\gamma y + \delta x y

という微分方程式系が得られる.

deSolve パッケージで解く

まずモデル式(ここでは "LVmodel" という名前にした)を定義する.

これを解くには ode という関数を使う.
ode の第一引数 y にはモデル式の初期値を与える. 引数 times にはシミュレーション時間を数列で与える. func は解きたい微分方程式("LVmodel"), parms はモデルのパラメータを与える.

library(deSolve)
LVmodel <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    
    dx   <- alpha*x - beta*x*y 
    dy    <- -gamma * y  + delta*x*y
    
    return(list(c(dx, dy)))
  })
}

pars  <- c(alpha=0.4,beta=0.6,gamma=0.3,delta=0.5)
yini  <- c(x = 100, y = 150)
times <- seq(0, 400, by = 1)
out   <- ode(yini, times, LVmodel, pars)

計算結果を図示する.

## ggplot2
library(ggplot2)
library(tidyr)
theme_set(theme_grey(20))
out2 <- gather(as.data.frame(out),variable, population, -time)
ggplot(out2,aes(x=time,y=population,colour=variable)) +
  geom_line(size=1.2)+labs(colour="")

f:id:abrahamcow:20150218142108p:plain

参考文献

微分方程式で数学モデルを作ろう

微分方程式で数学モデルを作ろう