commutateurs binaires dans desolve

Question

J'ai rencontré un problème en essayant de résoudre un ODE en r. J'ai un paramètre Q qui quand il devient plus grand que l'autre paramètre h arrête de couler dedans. l'ODE fonctionne très bien jusqu'à ce qu'il arrive au point où le commutateur se produit, puis arrête la course et me donne le 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 

le code est ci-dessous

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 
) 

Answer 1

Essayez d'utiliser une transition en douceur pour 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 définit la façon dont la frontière est forte.

Ainsi, votre modèle devrait ressembler à ceci:

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