Search code examples
rnandifferential-equations

NaNs produced in deSolve package

I have got a system of 8 differential equations that I am trying to solve using deSolve in R. It just returns NaN after the first few steps and doesn't solve it further. I tried various differential solvers like lsoda (default), bdf, adams, rk4 etc, but it didn't help.

Here is the sample R code:

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

I would like the X's be solved at least until they are reduced to 1% of the initial value. But, I see NaN values too early in the simulation. I tried changing the time step size to as low as 0.0005 and as high as 0.5 hr, but I still encounter the same problem.

Solution

  • These problems can be tough to diagnose, but I've started by refactoring your function - that is, I simplified the guts and made sure that it gave sufficiently close answers (using all.equal()) to the original.

    • your first 7 gradient values use identical forms, and can be collapsed into a single vectorized call. This is more efficient, easier to read, and easier to debug.
    • to make it easier to see where the problem is occurring, I further broke the expression into a series of simpler expressions (some of which are repeated in the original call: see #1 for advantages of reducing duplication ...)
    • factored out an unnecessary tmp3^2/tmp3 (== tmp3) in the equation
    • I put an if(any(is.na(grad))) test and a browser() call into the function so we can stop when the first NA/NaN value occurs and see what's going on ...
    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

    Now try running

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

    This drops us into an interactive browser environment.

    First intermediate expression looks OK (large, but finite):

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

    The second expression (3*state[v]/tmp1) looks OK, but note the last value is negative - this is presumably because the last (seventh) state variable has gone slightly negative.

     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

    Now when we try to take the cube root things go bad: unless a value is defined as a complex type, fractional powers of negative numbers are NaN in R

    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 will quickly spread through and mess up the entire state.

    At this point we could try track backward a bit farther and try to understand the floating-point inaccuracy that led to the negative value in the first place; however, it may or may not be easy (or even possible) to rewrite the expressions in a way that makes them sufficiently stable. This question and this question discuss ways of constraining solutions of ODEs to non-negative values ... the simplest (if it makes sense for your problem) is to put in a pmax(tmp3,0) or pmax(tmp3,very_small_positive_number) call to prevent negative values ...

    Minor comments:

    • are the factors of 3.14 in your questions supposed to be pi? R has a built-in pi value ...
    • for a given set of parameters, it looks like tmp1 is constant over time. You might want to pre-compute it outside of the gradient function for efficiency ...

    To see what was going on I added na.rm=TRUE to the sum as you suggested. I switched to method="euler"; this is less efficient, but simpler since it does very few intermediate calculations between gradient computations.

    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()

    enter image description here

    This shows that one component after another is dropping to very small values (and becoming NaN after we try to take the cube root of a negative value ...) Depending on what you were doing, it might be safe to set gradients whose values were NaN to zero ...