deSolve 包中产生的 NaN

Question

我有一个包含 8 个微分方程的系统，我正在尝试使用 R 中的 deSolve 来求解。它只是 returns 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

Problem/Question

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

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 and 讨论将 ODE 的解限制为非负值的方法...最简单的（如果它对您的问题有意义）是放入 pmax(tmp3,0) 或 pmax(tmp3,very_small_positive_number) 调用以防止负值 ...

次要评论：

你问题中3.14的因数应该是圆周率吗？ 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 的梯度设置为零是安全的...

deSolve 包中产生的 NaN

NaNs produced in deSolve package

r

nan

differential-equations

Problem/Question