使用 deSolve-Package 求解欧拉积分答案

【问题标题】：Using deSolve-Package to solve Euler integration使用 deSolve-Package 求解欧拉积分
【发布时间】：2016-09-12 17:54:27
【问题描述】：

我正在阅读 Soetaert 和 Herman 的“生态建模实用指南”。我以示例“6.7.2 Nutrient-Algae Chemostat Model - Euler Integration”的示例到达第 6 章 (Download)。因此，我尝试使用 deSolve-package 和 documentation 使用 R 复制建模。我完成了第一部分并得出了以下代码：

require(deSolve)

# Setting parameter values
parameters <- c(Nin=10, pmax=1, ks=1, dil=0.24)

# Setting state variable values
state <- c(DIN=0.1, Phyto=1)

# Defining function
Lorenz<-function(t, state, parameters) {
  with(as.list(c(state, parameters)),{
    # Rate of change
    dDIN <- -pmax * (DIN/DIN+ks) * Phyto - dil * (DIN-Nin)
    dPhyto <- pmax * (DIN/DIN+ks) * Phyto - dil * Phyto

    # Return the rate of change
    return(list(c(dDIN, dPhyto)))
  })
}

# Define time frame
times<-seq(0,20,by=0.1)

# Solving equations
out<-ode(y=state,times=times,func=Lorenz,parms=parameters,method="euler")
head(out)

# Plotting results
par(mfrow=c(1,3), oma = c(0, 0, 3, 0))
plot(out, xlab = "time", ylab = "-")
plot(out[, "Phyto"], out[, "DIN"], pch = ".")
mtext(outer = TRUE, side = 3, "Lorenz model", cex = 1.5)

所以我没有得到任何错误，但我知道这个解决方案是错误的。当我使用文档的示例代码时，一切正常。我的想法是错误在定义功能部分内。不将时间步长操作包含为 ...[i+1]... 似乎很奇怪，但是没有它们，文档也可以工作。

有没有人对deSolve 或整个练习有任何经验，可以给我任何解决问题的提示吗？

【问题讨论】：

标签： r ode

【解决方案1】：

你的错误是错误的括号，你必须使用这个：

dDIN <- -pmax * DIN/(DIN+ks) * Phyto - dil * (DIN-Nin)
dPhyto <- pmax * DIN/(DIN+ks) * Phyto - dil * Phyto

我想很清楚为什么，不是吗？

【讨论】：