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