1
$\begingroup$

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?

$\endgroup$
2
$\begingroup$

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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ @MichaelE2 As the resident expert on NDSolve, what are your thoughts? Thanks. $\endgroup$ – bbgodfrey Jun 8 at 19:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.