Search code examples
rif-statementdifferential-equationstimedelaydesolve

Error delay differential equation deSolve (dede)


I am writing a delay differential equation in deSolve (R), and I am getting an error message I am unsure how to solve. So for some background. I have a system with 12 differential equations, and 3 of those have a delay. I was successful in writing the system without deSolve, but I want to use deSolve because it gives me opportunity to easily work with other methods that are not a fixed-step euler. However, now I put it in deSolve, and I am using the general solver for delay differential equations (dede), I get an error.

This is the delay and the differential equation that the error is about:

lag2=ifelse(t-tau2<0, 0, e2(lagvalue(t-tau2,3))*lagvalue(t-tau2,8))

dn9dt=lag2-IP2*ethaP2M*P2-mu2*sf0B(B2)*P2-DNB(N2,B2)*P2+DD*PD

The first and third delay differential equations are done the same as this one and don't seem to have an error. The error is:

Error in lagvalue(t - tau2, 3) : illegal input in lagvalue - lag, 0, too large, at time = 15.945

Important to note is that the delay (tau2) is 16 in this case, and the error occurs sligtly before time = 16.

I have already tried to change t-tau2<0 to t-tau2<=0, but this doesn't help, and I have tried to increase the size of the history array (control=list(mxhist = 1e6)), which also didn't help. I have also tried to rewrite the lag a couple of times, but everytime I get the same error.

I have tried searching online, but I can hardly find anything on dede in deSolve, so I hope someone here might be able to help.


Solution

  • The question did not contain a full reproducible example, but the error can be reproduced if a simulation is run with a large number of time steps and the history array is too small. The following example is an adapted version from the ?dede help page:

    library("deSolve")
    
    derivs <- function(t, y, parms) {
    
      #cat(t, "\n") # uncomment this to see when the error occurs
      
      lag1 <- ifelse(t - tau1 < 0, 0.1, lagvalue(t - tau1, 2))
      lag2 <- ifelse(t - tau2 < 0, 0.1, lagvalue(t - tau2, 2))
    
      dy1 <- -y[1] * lag1 + lag2
      dy2 <-  y[1] * lag1 - y[2]
      dy3 <-  y[2] - lag2
      list(c(dy1, dy2, dy3))
    }
    
    yinit <- c(x=5, y=0.1, z=1)
    times <- seq(0, 40, by = 0.1)
    
    tau1 <-  1
    tau2 <- 10
    

    It runs successfuly with:

    yout <- dede(y = yinit, times = times, func = derivs, parms = NULL)
    

    but if we increase the number of time steps:

    times <- seq(0, 40, by = 1e-3)
    yout <- dede(y = yinit, times = times, func = derivs, parms = NULL)
    

    we can get an error:

    Error in lagvalue(t - tau2, 2) : 
      illegal input in lagvalue - lag, 0, too large, at time = 9.99986 
    

    It occurs after exceeding the threshold time (uncomment the cat above) when the algorithm starts to interpolate in the history array. As a solution, increase the history buffer:

    yout <- dede(y = yinit, times = times, func = derivs, 
      parms = NULL, control = list(mxhist = 1e5))
    

    It may also help to reduce the number of time steps.