binary switches in desolve

Question

I have run up against a problem in trying to solve a ODE in r. I have a parameter Q that when it gets bigger than the another parameter h stops having flow into it. the ODE works fine until it gets to the point where the switch occurs and then stops running and gives me the message:

  DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
  taken on this call before reaching TOUT     
  In above message, I1 = 5000

  In above message, R1 = 6.65299

 Warning messages:
  1: In lsoda(y, times, func, parms, ...) :
    an excessive amount of work (> maxsteps ) was done, but integration was 
   not successful - increase maxsteps
  2: In lsoda(y, times, func, parms, ...) :
    Returning early. Results are accurate, as far as they go

the code is below

parameters <- c(
                a = 0.032,
                b = (9 / 140),
                c = (5 / 1400),
                d = (95 / 700),
                k = 1 / 140,
                i = 0.25,
                # r = 0.2,
                n = 6000000,
                x = 0.5 ,
                y = 0.4,
                t = 1 / 180,        # important in looking at the shape
                u = 1 / 180,        # important in looking at the shape
                v = 1 / 360,        # important in looking at the shape
                p = 10,
                s = 10000,
                g = 100

                # e = .4,
                #h = 1000
)


state <- c(
            S = 5989900,
            E = 0,
            I = 0,
            Q = 0,
            D = 100,
            B = 0,
            C = 100,
            Y = 100,
            H = 1000,
            R = 1000,
            J = 1000,
            h = 1000,
            e = 0.1,
            r = 0.1
  )



 equation <- (function(t, state, parameters)
  with(as.list(c(state, parameters)), {
    # rate of change

    dS <- (-(a * S * I) / n) - (((1 / r) * S * D) / n)
    dE <- (a * S * I) / n + (((1 / r) * S * D) / n) - i * E
    if (h > Q)
      j = 1
    else if (h <= Q) 
      j = 0
    dI <- i * (j) * E - (e) * I - c * I - d * I
    dQ <- (j) * (e) * I - b * Q - k * Q
    dD <- d * I - r * D
    dB <- b * Q + r * D
    dC <- c * I + k * Q

    dY <- p * (b * Q + r * D)
    dR <- (1 - x - y) * (p * (b * Q + r * D))  - t * (R)
    de <- t * (s / R)
    dJ <- (y) * (p * (b * Q + r * D))  - v * (J)
    dr <- v * (s / J)
    dH <- (x) * (p * (b * Q + r * D)) - u * (H)
    dh <- u * (H / g)


    # return the rate of change
    list(c(dS, dE, dI, dQ, dD, dB, dC, dY, dR, de, dJ, dr, dH, dh))
  }))
#

# solve the equations for certain starting parameters


library(deSolve)
times <- seq(0, 200, by = 1)

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

First, try increasing the maximum number of steps. `out <- ode(y = state, times = times, func = equation, parms = parameters, maxsteps = 1e5)`. If that doesn't work, think about using a function that smooths the transition somewhat instead of a having a hard boundary, which numerical solvers can have a problem with. — Dan, May 17 '17 at 11:09
Lyngbakr how would you go about doing smoothing the transition? Do you have any examples? I just tried increasing the max steps and it didn't work unfortunately. — angusr, May 17 '17 at 11:25
Ahhhh thanks Lyngbakr, I have been trying to plum it into the ODE, I keep getting the message: Error in checkFunc(Func2, times, y, rho) : The number of derivatives returned by func() (10012) must equal the length of the initial conditions vector (14). but I still have the same amount of derivatives? — angusr, May 18 '17 at 09:10
I've updated my answer with what I think your model should look like. Does this work? — Dan, May 18 '17 at 11:06
In an answer to [this question](https://stackoverflow.com/questions/61484483/r-ode-function-desolve-package-change-the-value-of-a-parameter-as-a-function?rq=1) it was suggested to have a look at [deSolve help](https://tpetzoldt.github.io/deSolve-forcing/deSolve-forcing.html) pages on `forcing` and `events` which I've found particularly useful — Alf Pascu, May 21 '22 at 00:34

Dan · Accepted Answer · 2017-05-18T11:06:00.227

Try using a smoother transition for j:

library(sigmoid)

# Function will transition between 0 and 1 when h and Q are approximately equal
smooth.transition <- function(h, Q, tune = 0.01){
  sigmoid((h/Q - 1)/tune)
}

Q <- 1
h <- seq(0.001, 5, by = 0.001)
j <- smooth.transition(h, Q)

plot(h/Q, j, type = "l")

tune defines how sharp the boundary is.

So, your model should look something like this:

 equation <- (function(t, state, parameters
    with(as.list(c(state, parameters)), {
    # rate of change

    dS <- (-(a * S * I) / n) - (((1 / r) * S * D) / n)
    dE <- (a * S * I) / n + (((1 / r) * S * D) / n) - i * E

    ############################
    j <- smooth.transition(h, Q)
    ############################

    dI <- i * (j) * E - (e) * I - c * I - d * I
    dQ <- (j) * (e) * I - b * Q - k * Q
    dD <- d * I - r * D
    dB <- b * Q + r * D
    dC <- c * I + k * Q

    dY <- p * (b * Q + r * D)
    dR <- (1 - x - y) * (p * (b * Q + r * D))  - t * (R)
    de <- t * (s / R)
    dJ <- (y) * (p * (b * Q + r * D))  - v * (J)
    dr <- v * (s / J)
    dH <- (x) * (p * (b * Q + r * D)) - u * (H)
    dh <- u * (H / g)


    # return the rate of change
    list(c(dS, dE, dI, dQ, dD, dB, dC, dY, dR, de, dJ, dr, dH, dh))
}))

binary switches in desolve

1 Answers1