在我的颂歌模型中变量总是等于 0

Question

我使用

ode()

包中的

deSolve

函数在 R 中编写了一个 ODE 模型。该方程预测了叶面农药质量的演变。

app1

等于施用的杀虫剂质量（在第 156 天），函数应跟踪毒物的演变。但是，我注意到变量

mf

（叶质量）的输出始终等于 0，即使它应该随时间变化（从第 156 天开始）。我不确定是什么导致了这个问题。有人可以帮助我了解我做错了什么吗？

library(deSolve)

mod <- function(t, yini, param, Precip) {   
  app1 <- ifelse(t == 156, 1.12, 0)
  jgrow <- 32
  igrow <- ifelse(t >= 136 & t < 169, t - 136, ifelse(t > 168, 0.9, ifelse(t < 136, 0, NA)))
  cover <- 0.9 * igrow / jgrow
  P = Precip[t]
  
  with(as.list(c(yini, param)), { 
    
    mharv = ifelse(t== 201, mf, 0)
    
    
    d_MF <- (app1 * sa * (1-drift)*cover*sa)- (mf* exp(-kf*dt)) + (yf * (mf* exp(-kf*dt))) - mf * (1-exp(-fet * P * dt)) - mharv
    return(list(c(d_MF))) 
    
  })} 

params <- list(drift = 0.1, dt = 1, kf = 0.023, sa = 32.4, fet = 0.2, yf = 0.6)
t <- seq(1, 365, by = 1)
Precip=rep(0, 365)
Precip[155:205] =c(0, 0.96, 0.71, 0, 0, 0.08, 14.99, 0.2, 0.2, 0.25, 0, 0, 0.1, 0, 3.66, 0, 0.13, 0, 5.28, 0.03, 0, 0, 3.15, 0.84, 1.83, 0.03, 0, 2.29, 0, 0, 0.15, 1.65, 0, 0.91, 0.18, 2.01, 0, 0, 0, 0, 0, 3.17, 0.1, 0, 1.45, 0.25, 0, 0.03, 0.46, 0.03, 0)
y0 <- c(mf = 0)
sol <- ode(y = y0, times = t, func = mod, parms = params, Precip = Precip)

Answer 1

app1

和

mharv

实际上总是零，斜率函数几乎肯定不会在这个特定点被评估，即使被评估，这个单点尖峰的影响也可以忽略不计。

出于同样的原因，

t

几乎肯定不会在对

mod

的调用中成为整数，当您尝试使用

t

作为整数索引时应该会出错，但也许它会被自动截断。

所以最后，导数表达式中的每一项都有一个为零的因子，使零函数成为正确的解决方案。

Answer 2

这里是解决方案的一部分。时间用

floor

四舍五入，

precip

转换为具有恒定插值的强制函数。

igrow <- ifelse(...

函数看起来很奇怪，尤其是模型中的值永远不应该是 NA。

library(deSolve)

mod <- function(t, yini, param, f_precip) {   
  app1 <- ifelse(floor(t) == 156, 1.12, 0)   # use round, floor or trunc
  jgrow <- 32
  
  ## the following can also be converted to a forcing function
  ## but I suspect it may have some errors
  igrow <- ifelse(t >= 136 & t < 169, t - 136, ifelse(t > 168, 0.9, ifelse(t < 136, 0, NA)))
  
  cover <- 0.9 * igrow / jgrow
  P <- f_precip(t)
  
  with(as.list(c(yini, param)), { 
    mharv <- ifelse(floor(t) == 201, mf, 0) # use floor, trunc or round
    d_MF <- (app1 * sa * (1 - drift) * cover * sa) - (mf * exp(-kf * dt)) + 
            (yf * (mf * exp(-kf * dt))) - mf * (1 - exp(-fet * P * dt)) - mharv
    return(list(c(d_MF))) 
  })
} 

params <- list(drift = 0.1, dt = 1, kf = 0.023, sa = 32.4, fet = 0.2, yf = 0.6)

t <- seq(1, 365, by = 1)
Precip <- rep(0, 365)
Precip[155:205] <- c(0, 0.96, 0.71, 0, 0, 0.08, 14.99, 0.2, 0.2, 0.25, 0, 0, 0.1, 
                     0, 3.66, 0, 0.13, 0, 5.28, 0.03, 0, 0, 3.15, 0.84, 1.83, 0.03, 
                     0, 2.29, 0, 0, 0.15, 1.65, 0, 0.91, 0.18, 2.01, 0, 0, 0, 0, 0, 
                     3.17, 0.1, 0, 1.45, 0.25, 0, 0.03, 0.46, 0.03, 0)

## define a forcing function by tabular interpolation
f_precip <- approxfun(t, Precip, method = "constant", rule = 2)

y0 <- c(mf = 0)
sol <- ode(y = y0, times = t, func = mod, parms = params, 
           f_precip = f_precip)

plot(sol)