星期四, 四月 14, 2011

用R语言来解微分方程模型

20世纪40年代,Lotka(1925)和Volterra(1926)奠定了种间竞争关系的理论基础,他们提出的种间竞争方程对现代生态学理论的发展有着重大影响。Lotka-Volterra模型由一个微分方程组来描述。下面我们就用R中的deSolve包来解这个方程组,代码如下:

#首先加载包
library(deSolve)

#设置模型参数
pars <- c(rI = 0.2, rG = 1.0, rM = 0.2, AE = 0.5, K = 10)
#设置模型的初值和时间
yini <- c(P = 1, C = 2)
times <- seq(0, 80, by = 1)
#构建微分方程,以作为求解的输入
LVmod0D <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
  IngestC <- rI * P * C
  GrowthP <- rG * P * (1 - P/K)
  MortC <- rM * C
  dP <- GrowthP - IngestC
  dC <- IngestC * AE - MortC
  return(list(c(dP, dC)))
  })
  }

#利用ode函数求解结果
out<- ode(func=LVmod0D,y=yini,parms=pars,times=times)
#利用绘图函数画出时序图
matplot(out[,"time"], out[,2:3], type = "l", xlab = "time", ylab = "Conc", main = "Lotka-Volterra", lwd = 2)
legend("topright", c("prey", "consumer"), col = 1:2, lty = 1:2)

没有评论:

发表评论