deSolve ODE 不适用于微分方程（计算 NA）
deSolve ODE Not Working with Differential Equations (Calculates NA)
我一直在使用方法 deSolve::ode45，该方法一直有效，直到我对方程式进行了一些必要的更改。有谁知道为什么 ODE 求解器不工作？我尝试了 运行 ode45 以及默认的 ode 方法，但均无效。如果有任何进一步的解释会有帮助，请告诉我。
我检查了微分方程，我相信它们是正确的。
使用的方程式如下：
CCHFModel = function(t,x,params)
{
# get SIR values
SH <- x[1]
EH <- x[2]
IA <- x[3]
IS <- x[4]
RH <- x[5]
ST <- x[6]
IT <- x[7]
SC <- x[9]
IC <- x[10]
RC <- x[11]
# Load values ----
# Beta values
betaHHA = params["betaHHA"]
betaHHS = params["betaHHS"]
betaTH = params["betaTH"]
betaCH = params["betaCH"]
betaTC = params["betaTC"]
betaCT = params["betaCT"]
betaTT = params["betaTT"]
# Gamma value
gamma = params["gamma"]
# death rates
muH = params["muH"]
muT = params["muT"]
muC = params["muC"]
# birth rates
piH = params["piH"]
piT = params["piT"]
piC = params["piC"]
# incubation
deltaHS = params["deltaHS"]
deltaHA = params["deltaHA"]
# recovery rate
alphaA = params["alphaA"]
alphaS = params["alphaS"]
alphaC = params["alphaC"]
# total population
NH = (SH + IA + IS + EH + RH) + (piH * SH) - (muH * SH)
NT = (ST + IT) + (piT * ST) - (muT * ST)
NC = (SC + IC + RC) + (piC * SC) - (muH * SC)
# tick carrying Capacity
# KT = NC * 130 # 130 ticks per carrier max
#computations ----
dSHdt <- (piH * NH) - (betaHHA * IA + betaHHS * IS + betaCH * IC + betaTH * IT)*(SH/NH) - (muH * SH)
dEHdt <- (betaHHA * IA + betaHHS * IS + betaCH * IC + betaTH * IT)*(SH/NH) - ((deltaHA + muH)*EH)
dIAdt <- (deltaHA * EH) - ((alphaA + muH + deltaHS) * IA)
dISdt <- (deltaHS * IA) - ((alphaS + muH + gamma) * IS)
dRHdt <- alphaA * IA + alphaS * IS - muH*RH
dSTdt <- (piT * NT) - (betaTT * IT + betaCT * IC)*(ST/NT) - (muT * ST)
dITdt <- (betaTT * IT + betaCT * IC)*(ST/NT) - (muT * IT)
dSCdt <- (piC * NC) - (betaTC * IT)*(SC/NC) - (muC * SC)
dICdt <- (betaTC * IT)*(SC/NC) - ((alphaC +muC) * IC)
dRCdt <- (alphaC * IC) - (muC * RC)
# return results
list(c(dSHdt, dEHdt, dIAdt, dISdt, dRHdt, dSTdt, dITdt, dSCdt, dICdt, dRCdt))
}
我 运行 ODE 求解器使用：
defaultParms = c(betaHHA = .0413,
betaHHS = .0413,
betaTH = .2891,
betaCH = .0826,
betaTC = (1/365),
betaCT = 59/365,
betaTT = ((1/(365 * 2)) * .04) * 280,
gamma = 1/10,
muH = (1/(365 * 73)),
muT = (1/(365 * 2)),
muC = (1/(11 * 365)),
piH = 1.25/(73 * 365),
piT = 4.5/730,
piC = 1/(11 * 365),
deltaHS = 1/3,
deltaHA = 1/2,
alphaA = 1/17,
alphaS = 1/17,
alphaC = 1/7)
# time to start solution
t = seq(from = 0, to = 365, by = 0.1)
#initialize initial conditions
initialConditions = c(SH = 10000, EH = 5, IA = 5, IS = 10, RH = 2, ST = 80000, IT = 50, SC = 30000, IC = 5, RC = 1)
dataSet = ode(y = initialConditions, times = t, func = CCHFModel, parms = defaultParms)%>%
as.data.frame()
在 运行 之后，所有遵循初始条件的输出都是 NA。
这是由于打字错误 - 您在代码的第一部分错误编号了输入值的翻译（即，您跳过了 x[8]。我将完成两个（希望如此）有用的练习，首先解释我是如何调试它然后展示如何重写你的函数以减少它的error-prone ...
调试
- 尝试 运行
t=0 的梯度函数，x=<initial conditions>:
CCHFModel(0,initialConditions, defaultParms)
## piH betaHHA deltaHA deltaHS alphaA piT
## -15.02882327 12.62349834 0.53902803 0.07805607 0.88227788 385.31052332
## betaTT piC betaTC alphaC
## 0.85526763 NA NA NA
嗯，我们已经看到我们遇到了问题。为什么计算梯度的最后三个元素是 NA?
在函数末尾附近（dsCdt <- ... 行之前）添加 browser() 以便我们仔细查看。重新定义函数，再次尝试计算梯度。
当我们到达那里并打印出计算中涉及的一些量时，我们看到 NC 和 RC 都是 NA ... 我们也可以看出 RC 的 NA 值会导致 NC 成为 NA，所以让我们检查一下 RC ...[=45= 的定义]
啊哈！ RC定义为x[11]，但是length(initialConditions)只有10……仔细一看，我们漏掉了x[8]。重新定义正确地给出了非NA的值（我不知道它们是否正确，但至少它们不是NA）。
error-proofing (1)
虽然使用 [] 或 [[]] 来提取向量的元素通常会给出相同的答案，当您想提取单个元素时，您应该始终使用 [[]]向量中的元素（标量）。原因如下：
initialConditions[11] ## NA
initialConditions[[11]] ## Error in x[[11]] : subscript out of bounds
如果您使用 []，NA 会通过您的代码传播，您必须寻找原始来源。如果您使用 [[]]，R 会立即失败并告诉您问题出在哪里。另一个好处是 [] 以一种通常没有意义的方式传播向量元素的名称（看看上面“debugging/1”中输出的名称......）
error-proofing (2)
您可以通过将解包代码（计算总人口之前的所有内容）替换为
来避免参数和状态向量的所有繁琐和 error-prone 解包
comb <- c(as.list(x), as.list(params))
attach(comb)
on.exit(detach(comb))
假设您的参数和状态向量已正确命名（并且它们之间没有重叠的名称），这将创建一个命名列表并允许在您的函数中按名称查找元素； on.exit(detach(comb)) 确保最后一切都得到正确清理。 （您会看到使用 with() 执行此操作的建议；我更喜欢此处的策略，因为它使函数内的调试 [如有必要] 更容易。但正如@tpetzoldt 在评论中指出的那样，您应该 始终 将 attach(...) 与 on.exit(detach(...)) 配对；否则事情会变得非常混乱和混乱...)
在函数的末尾我会使用
g <- c(dSHdt, dEHdt, dIAdt, dISdt, dRHdt, dSTdt, dITdt, dSCdt, dICdt, dRCdt)
names(g) <- names(x)
list(g)
确保梯度向量被正确标记，这使得故障排除更容易。
我一直在使用方法 deSolve::ode45，该方法一直有效，直到我对方程式进行了一些必要的更改。有谁知道为什么 ODE 求解器不工作？我尝试了 运行 ode45 以及默认的 ode 方法，但均无效。如果有任何进一步的解释会有帮助，请告诉我。
我检查了微分方程，我相信它们是正确的。
使用的方程式如下：
CCHFModel = function(t,x,params)
{
# get SIR values
SH <- x[1]
EH <- x[2]
IA <- x[3]
IS <- x[4]
RH <- x[5]
ST <- x[6]
IT <- x[7]
SC <- x[9]
IC <- x[10]
RC <- x[11]
# Load values ----
# Beta values
betaHHA = params["betaHHA"]
betaHHS = params["betaHHS"]
betaTH = params["betaTH"]
betaCH = params["betaCH"]
betaTC = params["betaTC"]
betaCT = params["betaCT"]
betaTT = params["betaTT"]
# Gamma value
gamma = params["gamma"]
# death rates
muH = params["muH"]
muT = params["muT"]
muC = params["muC"]
# birth rates
piH = params["piH"]
piT = params["piT"]
piC = params["piC"]
# incubation
deltaHS = params["deltaHS"]
deltaHA = params["deltaHA"]
# recovery rate
alphaA = params["alphaA"]
alphaS = params["alphaS"]
alphaC = params["alphaC"]
# total population
NH = (SH + IA + IS + EH + RH) + (piH * SH) - (muH * SH)
NT = (ST + IT) + (piT * ST) - (muT * ST)
NC = (SC + IC + RC) + (piC * SC) - (muH * SC)
# tick carrying Capacity
# KT = NC * 130 # 130 ticks per carrier max
#computations ----
dSHdt <- (piH * NH) - (betaHHA * IA + betaHHS * IS + betaCH * IC + betaTH * IT)*(SH/NH) - (muH * SH)
dEHdt <- (betaHHA * IA + betaHHS * IS + betaCH * IC + betaTH * IT)*(SH/NH) - ((deltaHA + muH)*EH)
dIAdt <- (deltaHA * EH) - ((alphaA + muH + deltaHS) * IA)
dISdt <- (deltaHS * IA) - ((alphaS + muH + gamma) * IS)
dRHdt <- alphaA * IA + alphaS * IS - muH*RH
dSTdt <- (piT * NT) - (betaTT * IT + betaCT * IC)*(ST/NT) - (muT * ST)
dITdt <- (betaTT * IT + betaCT * IC)*(ST/NT) - (muT * IT)
dSCdt <- (piC * NC) - (betaTC * IT)*(SC/NC) - (muC * SC)
dICdt <- (betaTC * IT)*(SC/NC) - ((alphaC +muC) * IC)
dRCdt <- (alphaC * IC) - (muC * RC)
# return results
list(c(dSHdt, dEHdt, dIAdt, dISdt, dRHdt, dSTdt, dITdt, dSCdt, dICdt, dRCdt))
}
我 运行 ODE 求解器使用：
defaultParms = c(betaHHA = .0413,
betaHHS = .0413,
betaTH = .2891,
betaCH = .0826,
betaTC = (1/365),
betaCT = 59/365,
betaTT = ((1/(365 * 2)) * .04) * 280,
gamma = 1/10,
muH = (1/(365 * 73)),
muT = (1/(365 * 2)),
muC = (1/(11 * 365)),
piH = 1.25/(73 * 365),
piT = 4.5/730,
piC = 1/(11 * 365),
deltaHS = 1/3,
deltaHA = 1/2,
alphaA = 1/17,
alphaS = 1/17,
alphaC = 1/7)
# time to start solution
t = seq(from = 0, to = 365, by = 0.1)
#initialize initial conditions
initialConditions = c(SH = 10000, EH = 5, IA = 5, IS = 10, RH = 2, ST = 80000, IT = 50, SC = 30000, IC = 5, RC = 1)
dataSet = ode(y = initialConditions, times = t, func = CCHFModel, parms = defaultParms)%>%
as.data.frame()
在 运行 之后，所有遵循初始条件的输出都是 NA。
这是由于打字错误 - 您在代码的第一部分错误编号了输入值的翻译（即，您跳过了 x[8]。我将完成两个（希望如此）有用的练习，首先解释我是如何调试它然后展示如何重写你的函数以减少它的error-prone ...
调试
- 尝试 运行
t=0的梯度函数，x=<initial conditions>:
CCHFModel(0,initialConditions, defaultParms)
## piH betaHHA deltaHA deltaHS alphaA piT
## -15.02882327 12.62349834 0.53902803 0.07805607 0.88227788 385.31052332
## betaTT piC betaTC alphaC
## 0.85526763 NA NA NA
嗯，我们已经看到我们遇到了问题。为什么计算梯度的最后三个元素是 NA?
在函数末尾附近（
dsCdt <- ...行之前）添加browser()以便我们仔细查看。重新定义函数，再次尝试计算梯度。当我们到达那里并打印出计算中涉及的一些量时，我们看到
NC和RC都是NA... 我们也可以看出RC的NA值会导致NC成为NA，所以让我们检查一下RC...[=45= 的定义]啊哈！
RC定义为x[11]，但是length(initialConditions)只有10……仔细一看，我们漏掉了x[8]。重新定义正确地给出了非NA的值（我不知道它们是否正确，但至少它们不是NA）。
error-proofing (1)
虽然使用 [] 或 [[]] 来提取向量的元素通常会给出相同的答案，当您想提取单个元素时，您应该始终使用 [[]]向量中的元素（标量）。原因如下：
initialConditions[11] ## NA
initialConditions[[11]] ## Error in x[[11]] : subscript out of bounds
如果您使用 []，NA 会通过您的代码传播，您必须寻找原始来源。如果您使用 [[]]，R 会立即失败并告诉您问题出在哪里。另一个好处是 [] 以一种通常没有意义的方式传播向量元素的名称（看看上面“debugging/1”中输出的名称......）
error-proofing (2)
您可以通过将解包代码（计算总人口之前的所有内容）替换为
来避免参数和状态向量的所有繁琐和 error-prone 解包comb <- c(as.list(x), as.list(params))
attach(comb)
on.exit(detach(comb))
假设您的参数和状态向量已正确命名（并且它们之间没有重叠的名称），这将创建一个命名列表并允许在您的函数中按名称查找元素； on.exit(detach(comb)) 确保最后一切都得到正确清理。 （您会看到使用 with() 执行此操作的建议；我更喜欢此处的策略，因为它使函数内的调试 [如有必要] 更容易。但正如@tpetzoldt 在评论中指出的那样，您应该 始终 将 attach(...) 与 on.exit(detach(...)) 配对；否则事情会变得非常混乱和混乱...)
在函数的末尾我会使用
g <- c(dSHdt, dEHdt, dIAdt, dISdt, dRHdt, dSTdt, dITdt, dSCdt, dICdt, dRCdt)
names(g) <- names(x)
list(g)
确保梯度向量被正确标记，这使得故障排除更容易。