NDSolveValue::ndinid: Initial condition Sign[x8\$...] is not on the range specified by the discrete variable NDSolves\$

Bug introduced in 12.3 or earlier and persisting through 12.3.1 or later

I've checked the question about NDSolve::ndinid

This is the equation(Sorry for the complexity of this equation) and it's already written by Piecewise and on the right side as Michael E2 mentioned. However, the problem still exist.

μm=0.67;
ks=0.28;
x2max=2039;
x3max=939.5;
x4max=1026;
x5max=360.9;
x7max=80;
kgm=0.53;
kpm=0.14;

NR[expr_]:=Piecewise[{{1,expr>0},{0,expr<= 0}}];
μ = (μm*x2[t])/(
x2[t] + ks) (1 - x7[t]/x7max) (1 - x2[t]/x2max) (1 -
x3[t - τ1]/x4max) (1 - x4[t]/x4max) (1 - x5[t]/x5max);
q20 = m2 + μ/y2 + dq2*x2[t]/(x2[t] + k2max);
q2 = (k1*x2[t])/(x2[t] + k2) +
k3 (x2[t] - x6[t - τ2])*NR[x2[t] - x6[t - τ2]];
q3 = l1*(k4*x8[t - τ4])/(x8[t - τ4] + k5) +
l2*k6*(x8[t - τ4] - x3[t - τ1])*
NR[x8[t - τ4] - x3[t - τ1]];
q4 = m4 + μ*y4;
q5 = m5 + μ*y5;
l1=1;
l2=1;

cond = {x1[t /; t <= 0] == 0.102, x2[t /; t <= 0] == 418.2609,
x3[t /; t <= 0] == 0, x4[t /; t <= 0] == 0, x5[t /; t <= 0] == 0,
x6[t /; t <= 0] == 0, x7[t /; t <= 0] == 0, x8[t /; t <= 0] == 0};
pars={k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,k17,m2,m4,m5,y2,y4,y5,dq2,k2max};
minPara =
v1 -> 74.15520144717327, v2 -> 20.601402223057743,
v3 -> 202.89918221914687, v4 -> 124.79211581419787,
v5 -> 2.4441045037299847, v6 -> 47.960088480721915,
v7 -> 9.272966935510649, v8 -> 55.015434647393015,
v9 -> 4.369360890981919, v10 -> 2103.346395453384,
v11 -> 8.387912952673362, v12 -> 268.3929837388525,
v13 -> 33.53996289971859, v14 -> -3.224365183029838,
v15 -> 11.193999240900595, v16 -> 33.05890696229195,
v17 -> 88.20748789196996, v18 -> -3.3751062489269783,
v19 -> -3.673156600729279, v20 -> 3.8923944544765714,
v21 -> -11.123011817112472, v22 -> 31.8854734957005,
v23 -> 11.36136104963341, v24 -> 55.43326765393106,
v25 -> 9.250141726334826}[[;; , 2]]];

eqn = {
x1'[t] == μ*x1[t],
x2'[t] == -q2*x1[t],
x3'[t] == q3*x1[t],
x4'[t] == q4*x1[t],
x5'[t] == q5*x1[t],
x6'[t] ==
1/k7*((k8*x2[t])/(x2[t] + k9) +
k10 (x2[t] - x6[t - τ2])*NR[x2[t] - x6[t - τ2]] -
q20) - μ*x6[t - τ2],
x7'[t] == (k11*x6[t - τ2])/(
kgm (1 + x7[t - τ3]/k12) + x6[t - τ2]) - (
k13*x7[t - τ3])/(
kpm + x7[t - τ3] (1 + x7[t - τ3]/k14)) - μ*
x7[t - τ3],
x8'[t] == (k13*x7[t - τ3])/(
kpm + x7[t - τ3] (1 + x7[t - τ3]/k14)) -
l1*k15*x8[t - τ4]/(x8[t - τ4] + k16) -
l2*k17 (x8[t - τ4] - x3[t - τ1]) NR[
x8[t - τ4] - x3[t - τ1]] - μ*x8[t - τ4]
}/. minPara;

tauNDSolve[v1_, v2_, v3_, v4_] :=
NDSolveValue[{eqn /. {τ1 -> v1, τ2 -> v2, τ3 ->
v3, τ4 -> v4}, cond}, {x1, x2, x3, x4, x5, x6, x7,
x8}, {t, 0, 7}];

The error is

NDSolveValue::ndinid: Initial condition Sign[x8\$...] is not on the range specified by the discrete variable NDSolves\$..

Is there any modification to make it work? Grateful for any suggestion!

• There is probably some constant missing in the first DE, x1'[t] == *x1[t]. Jul 30 '21 at 16:47
• Sorry, it's μ * x1[t]. Jul 30 '21 at 17:49
• I think you should report it to WRI. I’m not sure “delayed” events (such as represented by your Piecewise) are supported. There’s nothing about it in the docs that I’ve found that indicate whether it is or isn’t. I was unable to get even a simple system with a delayed event to work. If unsupported, better parsing might give an error message to that effect. Jul 30 '21 at 23:36

The problem lies in your definition of NR. Apparently, Mathematica has some problems with these piecewise-like functions (see here and here).

Try changing it to:

NR[expr_ /; expr > 0] := 1;
NR[expr_ /; expr <= 0] := 0;

This works for me in Mathematica 12.3 - I can evaluate it at say $$(1,1,1,1$$).

tauNDSolve[1, 1, 1, 1] • Oh, right, I didn't notice it stopped at 2.4 ... Jul 30 '21 at 18:27
• That the computation fails at about 2.4 has nothing to do with the original issue raised in this question. Instead, it fails because 0.14 + (1 - 0.3101385678220016 x7[t - v3]) x7[t - v3] vanishes there, causing x7'[t] and x8[t]` to become singular. This is an inherent problem with the ODEs for the parameters chosen. Therefore, the answer above correctly identifies problem and provides an effective workaround. Jul 31 '21 at 15:26

As @bbgodfrey pointed out, the DDE in tauNDSolve[1, 1, 1, 1] has a singularity at t = ~2.4. This is coincidentally near a discontinuity event due to NR[], so I mistakenly thought that was the problem. But with that cleared up, my original idea of turning off discontinuity processing works:

tauNDSolve[v1_, v2_, v3_, v4_] :=
NDSolveValue[{eqn /. {\[Tau]1 -> v1, \[Tau]2 -> v2, \[Tau]3 ->
v3, \[Tau]4 -> v4}, cond}, {x1, x2, x3, x4, x5, x6, x7, x8}, {t,
0, 7}, Method -> "DiscontinuityProcessing" -> False]

sol = tauNDSolve[1, 1, 1, 1] I still think the error message, which contains internal local symbols, should be reported to WRI. That shouldn't happen.