ParametricNDSolve “delayed time” error message

Question

Issue reported to Wolfram, Inc as a possible bug in Version 12.1; CASE:4554034.

When trying to solve a system of delay-differential equations with ParametricNDSolve, Mathematica throws the error

"Computed delayed time -d$18077+t$18078 = 6.` is in advance of the \
current integration time t$18078 = 5"

followed by several other error messages. The code I use is

ClearAll[Sus, Inf, beta, d,  inf0, tmax, 
sol, Sus0, Inf0]
tmax = 10^4;(*tmax is the horizon of solving.*)
sol = 
 ParametricNDSolve[{Sus'[t] == 
    Piecewise[{{-beta*Sus[t]*Inf[t], t >= 0}, {0, t < 0}}], 
   Inf'[t] == 
    Piecewise[{{beta*Sus[t]*Inf[t], 
       0 <= t < d}, {beta*Sus[t]*Inf[t] - beta*Sus[t - d]*Inf[t - d], 
       t >= d}, {0, 0 > t}}], Sus[t /; t <= 0] == 1 - inf0, 
   Inf[t /; t <= 0] == inf0}, {Sus, Inf}, {t, -d, 
      tmax - d}, {{inf0, 0, 0.5}, {beta, 0.1, 1}, {d, 1, 10}}];
Sus0 = Sus /. First@sol;
Inf0 = Inf /. First@sol;
Table[Sus0[0.05, 0.6, 5][i], {i, 1, 20}]

The answers from Table[Sus0[...]] are wrong (negative).

I did not find anything about the error message on https://reference.wolfram.com/ or by web search.

How to obtain correct solutions from ParametricNDSolve, or failing that, more information about the delay time error?

Answer 1

Apparently, ParametricNDSolve does not correctly substitute d into Piecewise and the limits of integration. As a result, the wrong branches of Piecewise sometimes are chosen. Instead, use

tmax = 10^4;    
sol = With[{d = 5}, 
    ParametricNDSolve[{Sus'[t] == Piecewise[{{-beta*Sus[t]*Inf[t], t >= 0}, {0, t < 0}}], 
    Inf'[t] == Piecewise[{{beta*Sus[t]*Inf[t], 0 <= t < d}, {beta*Sus[t]*Inf[t] - 
        beta*Sus[t - d]*Inf[t - d], t >= d}, {0, 0 > t}}], 
    Sus[t /; t <= 0] == 1 - inf0, Inf[t /; t <= 0] == inf0}, 
    {Sus, Inf}, {t, -d, tmax - d}, {{inf0, 0, 0.5}, {beta, 0.1, 1}}]];
Sus0 = Sus /. First@sol;
Inf0 = Inf /. First@sol;
Table[Sus0[0.05, 0.6][i], {i, 1, 20}]

(* {0.912491, 0.85125, 0.758494, 0.632844, 0.486113, 0.345541, 0.234354, 
    0.157287, 0.108609, 0.0796167, 0.0628735, 0.0531666, 0.0473503, 0.0436512, 
    0.0410966, 0.0391615, 0.0375677, 0.0361705, 0.0348947, 0.033701} *)

Perhaps, this should be viewed as a bug. In any case, the documentation should be improved.