DDE example not going to zero as time -> large

[slwu89]

@sdwfrost in looking at the SDDE model issue https://github.com/SciML/StochasticDelayDiffEq.jl/issues/48 I decided to set the diffusion matrix to all zeros to see if the I state variable would decay to zero smoothly like a DDE. It does not; in fact this behavior can also be seen in the regular DDE model too, see example below, just copy-pasting the DDE tutorial in this repo, but solving to tmax = 1e4.

Checking the last state returned by the solver gives:

julia> sol_dde[end]
3-element Vector{Float64}:
   1.8036818718825726e-7
  10.002420763561224
 989.9975790560709

example:

using DifferentialEquations
using DelayDiffEq

function sir_dde!(du,u,h,p,t)
    (S,I,R) = u
    (β,c,τ) = p
    N = S+I+R
    infection = β*c*I/N*S
    (Sd,Id,Rd) = h(p, t-τ) # Time delayed variables
    Nd = Sd+Id+Rd
    recovery = β*c*Id/Nd*Sd
    @inbounds begin
        du[1] = -infection
        du[2] = infection - recovery
        du[3] = recovery
    end
    nothing
end;

δt = 0.1
tmax = 10000.0
tspan = (0.0,tmax)
t = 0.0:δt:tmax;

u0 = [990.0,10.0,0.0]; # S,I.R

function sir_history(p, t)
    [1000.0, 0.0, 0.0]
end;

p = [0.05,10.0,4.0]; # β,c,τ

prob_dde = DDEProblem(DDEFunction(sir_dde!),
        u0,
        sir_history,
        tspan,
        p;
        constant_lags = [p[3]]);

alg = MethodOfSteps(Tsit5());

sol_dde = solve(prob_dde,alg);

sol_dde[end]

[sdwfrost]

It's no coincidence here that I(t)->I(0) as t->infinity. Similar to your jump system, there appears to be a problem with the initial conditions to deal with the recovery of the 10 individuals that are all infected at t=0. I thought that one wouldn't have to deal with this specifically.

[slwu89]

I made an issue on the repo: https://github.com/SciML/DelayDiffEq.jl/issues/225

[devmotion]

I just checked quickly and observed the same behaviour with Mathematica. I'm not completely familiar with its NDSolve though, so I am not sure how well it deals with discontinuous history functions. Nevertheless, it makes me wonder if this delay SIR formulation should actually yield I(t) -> 0 as t -> Inf? More concretely, is I'(t) < 0 if I(t) = 10 for some t > 0? Otherwise, the dynamics would be limited to values I(t) >= 10 by the model.

[devmotion]

(I'm also a bit surprised that the recovery rate depends on the number of susceptibles in this model)

[sdwfrost]

To deal with the initial infecteds, one either needs to have a callback, or run an SI ODE model until t=tau and use that to initialize the DDE model at t=tau. I'll push an update shortly.

[sdwfrost]

Fixed (I think!) with https://github.com/epirecipes/sir-julia/commit/8f11553560e213734e9ee329a76d772810219123