deSolve封装NaNs

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小时，但我仍然遇到同样的问题。

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讨论了将ODEs的解约束为非负值的方法.最简单的(如果对您的问题有意义的话)是输入一个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的渐变设置为零可能是安全的。