在 deSolve 包中产生的 NaN

Question

我有一个由 8 个微分方程组成的系统，我试图在 R 中使用 deSolve 求解。 它只是在前几步之后返回NaN，并没有进一步解决它。 我尝试了各种微分求解器，例如lsoda（默认）、bdf、adams、rk4 等，但都没有帮助。

这里是示例 R 代码：

library(deSolve)
daero = c(5.29,4.16,2.49,1.53,0.7,0.41,0.21)*10^-4
rho = rep(1.27,7)
dgeo = daero * sqrt(1/rho)
r0 = dgeo/2
Fr = c(0.188,0.297,0.274,0.181,0.032,0.013,0.015)
X0 = Fr*200*10^-6
N0 = X0*(3/(4*3.14*r0^3*rho))

func1 <- function(t,state,parameters){
with(as.list(c(state,parameters)),{
dX1 = -D/((3*X1/(4*3.14*rho[1]*N0[1]))^(1/3))*(N0[1]*4*3.14* 
((3*X1/(4*3.14*rho[1]*N0[1]))^(1/3))^2)*(Cs-S/V)
dX2 = -D/((3*X2/(4*3.14*rho[2]*N0[2]))^(1/3))*(N0[2]*4*3.14* 
((3*X2/(4*3.14*rho[2]*N0[2]))^(1/3))^2)*(Cs-S/V)
dX3 = -D/((3*X3/(4*3.14*rho[3]*N0[3]))^(1/3))*(N0[3]*4*3.14* 
((3*X3/(4*3.14*rho[3]*N0[3]))^(1/3))^2)*(Cs-S/V)
dX4 = -D/((3*X4/(4*3.14*rho[4]*N0[4]))^(1/3))*(N0[4]*4*3.14* 
((3*X4/(4*3.14*rho[4]*N0[4]))^(1/3))^2)*(Cs-S/V)
dX5 = -D/((3*X5/(4*3.14*rho[5]*N0[5]))^(1/3))*(N0[5]*4*3.14* 
((3*X5/(4*3.14*rho[5]*N0[5]))^(1/3))^2)*(Cs-S/V)
dX6 = -D/((3*X6/(4*3.14*rho[6]*N0[6]))^(1/3))*(N0[6]*4*3.14* 
((3*X6/(4*3.14*rho[6]*N0[6]))^(1/3))^2)*(Cs-S/V)
dX7 = -D/((3*X7/(4*3.14*rho[7]*N0[7]))^(1/3))*(N0[7]*4*3.14* 
((3*X7/(4*3.14*rho[7]*N0[7]))^(1/3))^2)*(Cs-S/V)
dS = -dX1-dX2-dX3-dX4-dX5-dX6-dX7
list(c(dX1,dX2,dX3,dX4,dX5,dX6,dX7,dS))
})
}

state <- c(X1=X0[1],X2=X0[2],X3=X0[3],X4=X0[4],X5=X0[5],
       X6=X0[6],X7=X0[7],S=0)
parameters <- c(D=6.19*60*10^-6,rho=rho,N=N0,Cs=17*10^-6,V=1000)
times <- seq(0,3,by=0.0005)

out <- ode(y = state, times = times,func = func1, parms = parameters)
Output:
> out[1:20,]
    time           X1           X2           X3           X4
  0.0000 3.760000e-05 5.940000e-05 5.480000e-05 3.620000e-05
  0.0005 3.759491e-05 5.938700e-05 5.476652e-05 3.614143e-05
  0.0010 3.758982e-05 5.937400e-05 5.473305e-05 3.608290e-05
  0.0015 3.758473e-05 5.936100e-05 5.469959e-05 3.602440e-05
  0.0020 3.757964e-05 5.934800e-05 5.466613e-05 3.596594e-05
  0.0025 3.757455e-05 5.933500e-05 5.463268e-05 3.590750e-05
  0.0030 3.756947e-05 5.932201e-05 5.459924e-05 3.584910e-05
  0.0035 3.756438e-05 5.930901e-05 5.456581e-05 3.579074e-05
  0.0040 3.755929e-05 5.929602e-05 5.453238e-05 3.573240e-05
  0.0045 3.755420e-05 5.928303e-05 5.449897e-05 3.567410e-05
  0.0050 3.754912e-05 5.927004e-05 5.446556e-05 3.561583e-05
  0.0055 3.754403e-05 5.925705e-05 5.443215e-05 3.555760e-05
  0.0060 3.753895e-05 5.924406e-05 5.439876e-05 3.549940e-05
  0.0065 3.753386e-05 5.923107e-05 5.436537e-05 3.544123e-05
  0.0070 3.752878e-05 5.921808e-05 5.433199e-05 3.538310e-05
  0.0075 3.752369e-05 5.920510e-05 5.429862e-05 3.532499e-05
  0.0080 3.751861e-05 5.919212e-05 5.426525e-05 3.526692e-05
  0.0085 3.751352e-05 5.917913e-05 5.423190e-05 3.520889e-05
  0.0090          NaN          NaN          NaN          NaN
  0.0095          NaN          NaN          NaN          NaN

问题/疑问

我希望至少解决 X 直到它们减少到初始值的 1%。但是，我在模拟中看到NaN 值太早了。我尝试将时间步长更改为低至 0.0005 和高至 0.5 hr，但我仍然遇到同样的问题。

Answer 1

这些问题可能很难诊断，但我已经开始重构你的函数 - 也就是说，我简化了胆量并确保它给出了足够接近的答案（使用 all.equal()）原作。

• 您的前 7 个梯度值使用相同的形式，并且可以折叠成一个矢量化调用。这样更高效、更易于阅读且更易于调试。
• 为了更容易看出问题出在哪里，我进一步将表达式分解为一系列更简单的表达式（其中一些在原始调用中重复：请参阅 #1 以了解减少重复的优势...）
• 在等式中排除了不必要的 tmp3^2/tmp3 (== tmp3)
• 我在函数中添加了一个if(any(is.na(grad))) 测试和一个browser() 调用，这样我们就可以在第一个NA/NaN 值出现时停下来看看发生了什么...

func2 <- function(t,state,parameters, debug=TRUE){
    n <- length(state)
    v <- 1:(n-1)
    grad <- rep(NA,n)
    tmp1 <- (4*3.14*rho[v]*N0[v])
    tmp2 <- 3*state[v]/tmp1
    tmp3 <- tmp2^(1/3)
    grad[v] <- with(as.list(c(state,parameters)),{
        -D*(N0[v]*4*3.14*tmp3)*(Cs-S/V)
    })
    grad[n] <- -sum(grad[v])
    if (debug && any(is.na(grad))) browser()
    return(list(grad))
}
## test near-equality
all.equal(func1(0,state, parameters),func2(0,state, parameters)) ## TRUE

现在尝试运行

out <- ode(y = state, times = times,func = func2, parms = parameters)

这会将我们带入交互式浏览器环境。

第一个中间表达式看起来不错（大，但有限）：

Browse[2]> tmp1
[1]     8724442    28341529   121926846   347177124   640918307  1295801866
[7] 11127053948

第二个表达式 (3*state[v]/tmp1) 看起来不错，但请注意最后一个值是负数 - 这可能是因为最后一个（第七个）状态变量稍微变负了。

 Browse[2]> tmp2
           X1            X2            X3            X4            X5 
 1.289771e-11  6.262837e-12  1.333549e-12  3.037421e-13  2.588684e-14 
           X6            X7 
 3.751315e-15 -4.992697e-18 

现在，当我们尝试取立方根时，事情就变糟了：除非一个值被定义为 complex 类型，否则负数的分数幂在 R 中是 NaN

Browse[2]> tmp3
          X1           X2           X3           X4           X5           X6 
2.345151e-04 1.843276e-04 1.100702e-04 6.722049e-05 2.958192e-05 1.553798e-05 
          X7 
         NaN 

这会迅速蔓延并扰乱整个州。

在这一点上，我们可以尝试向后追溯一点，并尝试了解导致负值的浮点不准确性；但是，以使表达式足够稳定的方式重写表达式可能容易也可能不容易（甚至可能）。 This question 和 this question 讨论将 ODE 的解限制为非负值的方法...最简单的（如果对您的问题有意义的话）是调用 pmax(tmp3,0) 或 pmax(tmp3,very_small_positive_number) 以防止出现负值...

未成年人：

• 你的问题中 3.14 的因数应该是 pi 吗？ R 有一个内置的pi 值...
• 对于一组给定的参数，看起来tmp1 随着时间的推移是恒定的。您可能希望在梯度函数之外预先计算它以提高效率...

---

要查看发生了什么，我按照您的建议将na.rm=TRUE 添加到总和中。我切换到method="euler"；这效率较低，但更简单，因为它在梯度计算之间进行的中间计算非常少。

out <- ode(y = state, times = times,func = func2, parms = parameters,
           debug=FALSE,method="euler")
out <- out[rowSums(is.na(out))<9,]
png("SO_ode.png")
par(las=1)
matplot(out[,-1],type="l",lty=1,log="xy",col=1:8,
        xlab="time",ylab="")
dev.off()

这表明一个接一个的组件正在下降到非常小的值（并在我们尝试取负值的立方根后变为NaN...）根据您的操作，它可能 将值为NaN 的梯度设置为零是安全的...