r语言 - 处理ode时初始条件向量的导数个数和长度不等

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模型中，dA-方程中使用了变量BCG。由于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或整个互联网上搜索:

另一个注释:初始状态的数量级相当大(1e10)，这会导致数值问题。原则上，可以为每个状态变量设置单独的公差。然而，重新确定问题的规模通常更容易。我通常建议将其重新调整到0.0001到10000之间，或者在0.001到1000之间更好。