使用常微分方程时初始条件向量的导数和长度不等

Question

我对 R 中的 ODE 函数有疑问。我一直收到错误“

Error in checkFunc(Func2, times, y, rho) :  The number of derivatives returned by func() (3006) must equal the length of the initial conditions vector (6)

”。该函数似乎正在创建 tspan 大小加上初始条件向量数量的导数。关于如何进行的任何建议？

library(deSolve)
tend <- 300
dt <- 1
tspan <- seq(0, tend, by = dt)

# Parameters
rHI <- 0.4 # PI producers
rHN <- 0.4 # X2 free riders
rLI <- rHI * 0.1 # Xr resistant
rLN <- rHI * 0.1
K <- 10e13 # carrying capacity of tumor

delta0 <- 2
rhoS <- 1
rhoI <- 1
lambdaA <- 0.5
lambdaS <- 0.5
beta <- 1
eta <- 1e3
lambdaAS <- 0.02
BCG_dose <- 4e6
t0 <- 0      
tnow <- times[1]


## Initial Conditions
TV <- 1.4 * 10e10
HI0 <- 0.5 * TV
HN0 <- 0.0 * TV
LI0 <- 0.25 * TV
LN0 <- 0.25 * TV
A0 <- 0
S0 <- 0


##Setting Up 
initcond <- c(HI0, HN0, LI0, LN0, A0, S0)
names(initcond) <- c("HI", "HN", "LI", "LN", "A", "S")
##Setting Up BCG dose
interval <- 7
BCG <- rep(0, length(tspan))
BCG[seq(interval, length(tspan), by = interval)] <- BCG_dose
for (i in 2:length(tspan)) {
  if (tspan[i] - tnow >= interval) {
    # update tnow and BCG0 for exponential decay
    BCG0 <- BCG[i-1]
    t0 <- tnow
    tnow <- tspan[i]
    BCG[i] <- BCG0 * exp(-0.3 * (tnow - t0))
  } else {
    BCG[i] <- BCG[i-1]
  }
}
bladder_ODE <- function(t, x, params, BCG) {
  HI <- x[1]
  HN <- x[2]
  LI <- x[3]
  LN <- x[4]
  A <- x[5]
  S <- x[6]
  

  
  BCG0 <- BCG[which.max(tspan <= t)]
  dBCG <- -BCG0 * exp(-0.3 * (t - t0))
  BCG[which.max(tspan <= t)] <- BCG0 + dBCG
  
  delta <- delta0 * max((A/(S+A) - 0.5), 0)
  
  dHI <- rHI * HI * (1 - (HI + HN + LI + LN) / K) - delta * HI
  
  dHN <- rHN * HN * (1 - (HI + HN + LI + LN) / K)
  
  dLI <- rLI * LI * (1 - (HI + HN + LI + LN) / K)
  
  dLN <- rLN * LN * (1 - (HI + HN + LI + LN) / K)
  
  dA <- rhoS * (HI + LI) + beta * BCG / (eta + BCG) - lambdaA * A
  
  dS <- rhoI * (HN + LN) - lambdaS * S - lambdaAS * A * S
  
  return(list(c(dHI, dHN, dLI, dLN, dA, dS)))
}

out <- ode(y = initcond, times = tspan, func = bladder_ODE)

我已经尝试更改定义 tspan 的位置，但是，我仍然遇到同样的问题。我还将迭代的大小从 0.1 更改为 1，但是，它仍然不起作用。

Answer 1

问题是，在模型

bladder_ODE

中，变量

BCG

用于

dA

-方程。由于

BCG

是一个具有 301 个值的向量，因此方程被向量化并且

dA

被扩展到 301 个元素。

代码还包含另外两个小问题：

• 在第5行定义了

  tspan

  ，但是下面几行，

  tnow <- times[1]

  出现了，那可能是

  tnow <- tspan[1]

• 函数

  bladder_ODE <- function(t, x, params, BCG)

  包含一个参数

  BCG

  ，它也需要在

  ode

  函数调用中使用，见下文。

通过这些更改，现在可以重现错误。我在代码中放了两行

cat

行来打印受影响变量的长度。

bladder_ODE <- function(t, x, params, BCG) {
  HI <- x[1]
  HN <- x[2]
  LI <- x[3]
  LN <- x[4]
  A <- x[5]
  S <- x[6]

  BCG0 <- BCG[which.max(tspan <= t)]
  dBCG <- -BCG0 * exp(-0.3 * (t - t0))
  BCG[which.max(tspan <= t)] <- BCG0 + dBCG
  
  cat(length(BCG), "\n") # <----------
  
  delta <- delta0 * max((A/(S+A) - 0.5), 0)
  dHI <- rHI * HI * (1 - (HI + HN + LI + LN) / K) - delta * HI
  dHN <- rHN * HN * (1 - (HI + HN + LI + LN) / K)
  dLI <- rLI * LI * (1 - (HI + HN + LI + LN) / K)
  dLN <- rLN * LN * (1 - (HI + HN + LI + LN) / K)
  dA <- rhoS * (HI + LI) + beta * BCG / (eta + BCG) - lambdaA * A
  dS <- rhoI * (HN + LN) - lambdaS * S - lambdaAS * A * S
  
  cat(length(dA), "\n") # <-----------
  
  return(list(c(dHI, dHN, dLI, dLN, dA, dS)))
}

## aded argument BCG
out <- ode(y = initcond, times = tspan, func = bladder_ODE, BCG=BCG)

然后返回：

301 
301 
Error in checkFunc(Func2, times, y, rho) : 
  The number of derivatives returned by func() (306) must equal the length of the initial conditions vector (6)

除此之外，完整的

BCG[which.max(....)]

看起来很不寻常。您可以在 StackOverflow 或整个互联网上搜索：desolve forcings

另一条评论：初始状态的数量级非常大（1e10），这会导致数值问题。原则上，可以为每个状态变量设置单独的公差。但是，重新缩放问题通常更容易。我通常建议将其重新缩放到 0.0001 到 10000 之间的范围，甚至更好地在 0.001 到 1000 之间。