使用 deSolve-Package 求解欧拉积分
Using deSolve-Package to solve Euler integration
我正在研究 Soetaert 和 Herman 的 "Practical Guide to Ecological Modeling"。我到达了第 6 章(Download) with the example "6.7.2 Numerical Solution of a Nutrient-Algae Chemostat Model - Euler Integration". So I tried to replicate the modelling with R using the deSolve-package and the documentation。我完成了第一部分并导出了以下代码:
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
或整个练习有任何经验,可以给我任何解决问题的提示?
你的错误是错误的设置括号,你必须使用这个:
dDIN <- -pmax * DIN/(DIN+ks) * Phyto - dil * (DIN-Nin)
dPhyto <- pmax * DIN/(DIN+ks) * Phyto - dil * Phyto
我想原因很清楚了,不是吗?
我正在研究 Soetaert 和 Herman 的 "Practical Guide to Ecological Modeling"。我到达了第 6 章(Download) with the example "6.7.2 Numerical Solution of a Nutrient-Algae Chemostat Model - Euler Integration". So I tried to replicate the modelling with R using the deSolve-package and the documentation。我完成了第一部分并导出了以下代码:
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
或整个练习有任何经验,可以给我任何解决问题的提示?
你的错误是错误的设置括号,你必须使用这个:
dDIN <- -pmax * DIN/(DIN+ks) * Phyto - dil * (DIN-Nin)
dPhyto <- pmax * DIN/(DIN+ks) * Phyto - dil * Phyto
我想原因很清楚了,不是吗?