# Why do I get this nonsensical result for NDSolve with a delay-differential-algebraic system?

I have the following delay-differential-algebraic system:

\begin{align} c(t) =& \kappa\min (1, \tfrac{\beta}{\alpha} a(t)), \tag{1a}\\ a(t) =& \left\{ \begin{array}{lcl} \frac{\alpha}{\beta} & \text{if} & s(t) > 0,\\ \min(\frac{\alpha}{\beta}, c(t-\Delta)) & \text{if} & s(t) = 0 \end{array}\right. ,\tag{1b}\\ \frac{\text{d} s}{\text{d} t}(t) =& c(t-\Delta) - a(t).\tag{1c} \end{align}

which I am trying to numerically estimate by the following code:

ClearAll[a, c, t, s, inic, inia]
ClearAll["Global'*"]
kappa = 2;
alpha = 1;
beta = 2;
delta = 1;
amax = 10;
T = 1.1;
inic[t_] := -t;
inia[t_] := inic[t - delta]

eq1 = c[t] == kappa*Min[1, (beta/alpha)*a[t]]
eq2 = a[t] ==
Piecewise[{{Min[alpha/beta, c[t - delta]], s[t] == 0}, {alpha/beta,
s[t] > 0}}]
eq3 = s'[t] == c[t - delta] - a[t]

initialConditions = {
a[t /; t <= 0] == inia[t],
c[t /; t <= 0] == inic[t],
s[t /; t <= 0] == 0}
sol = NDSolve[
Union[{eq1, eq2, eq3}, initialConditions ], {a, c, s}, {t, -delta,
T}]

Plot[Evaluate[{c[t], a[t], s[t]} /. sol], {t, -delta, T}]


I don't know what I am doing wrong, but the plot does not make any sense: I've only included a small segment after $t=1$, because all three functions completely explode to $10^{20}$ and upwards at $t=2$. This is inconsistent with the equations, which state that $c(t), a(t)$ can attain at most values of $2$ and $\frac{1}{2}$ respectively, which they attain in the plot between $t=0$ and $t=1$. Moreover the plot doesn't show the values of the functions at interval $t=[-1,0]$, even though these are given as initial conditions.

Why does it give this strange result?

A piece of the puzzle may be that NDSolve gives the following error:

NDSolve::ivcon: The given initial conditions were not consistent with the differential-algebraic equations. NDSolve will attempt to correct the values. >>

• It appears the green line reaches 0 after t=1, causing a stiff system. – Feyre Oct 17 '16 at 10:00

I get two error messages, NDSolve::ivcon and NDSolve::ndsz, and a very strange (buggy?) set of interpolating functions.

Re NDSolve::ivcon:

subIC =   (* turn initial conditions into substitution rules *)
initialConditions /. {HoldPattern[Condition][p_, c_] :>
Inactive[Condition][Inactive[Pattern][p, Blank[]], c]} /.
{Equal -> Rule} // Activate
(*  {a[t_ /; t <= 0] -> 1 - t, c[t_ /; t <= 0] -> -t,  s[t_ /; t <= 0] -> 0}  *)

{eq1, eq2, eq3} /. t -> 0 /. subIC
(*  {False, False, Derivative[s] == 0}  *)


They're inconsistent because the two sides of eq1 and eq2 evaluate to different numbers:

List @@@ {eq1, eq2, eq3} /. t -> 0 /. subIC
(*  {{0, 2}, {1, 1/2}, {Derivative[s], 0}}  *)


Re NDSolve::ndsz: If you pay attention to the domain, you will see that interpolating functions have a domain of {-1, 1}. Plotting outside this domain (T = 1.1) leads to extrapolation.

Re strange functions:

They do not evaluate for t <= 0 (abnormal?):

{c[-0.5], a[-0.5], s[-0.5]} /. First@sol They evaluate between 0 and 0.9999 (normal):

Block[{t = {0.1, 0.5, 0.9}},
{c[t], a[t], s[t]} /. First@sol
]
(*
{{2., 2., 2.},
{0.5, 0.5, 0.5},
{0.045, 0.125, 0.045}}
*)


The extrapolation oscillates rapidly (normal, maybe?):

Block[{t = {1., 1.05, 1.08, 1.1}},
{c[t], a[t], s[t]} /. First@sol
]
(*
{{0., 1.18059*10^21, -4.72237*10^21, 0.},
{0., 2.95148*10^20, -1.18059*10^21, 0.},
{40.0292, 8.04778*10^19, 3.29637*10^20, 6.43821*10^20}}
*)


Replacement t -> val does not work and results in an invalid InterpolatingFunction[] (looks like a bug):

{c[t], a[t], s[t]} /. First@sol /. t -> 0.5


Function::flpar: Parameter specification {0.5} in Function[{0.5},{1-0.5,-0.5,0}] should be a symbol or a list of symbols.

Update: It appears that the NDSolveInitialHistory segment of the interpolating data has an expression that depends on the Global variable t:

(c[t] /. First@sol /. t -> 0.5)[[0, 4, 1]]


Function::flpar: Parameter specification {0.5} in Function[{0.5},{1-0.5,-0.5,0}] should be a symbol or a list of symbols.

(* note that the argument of  Function  has been replaced with 0.5
{2.44141*10^-8,
ExperimentalNumericalFunction[Function[{0.5}, {1 - 0.5, -0.5, 0}],
"-NumericalFunctionData-"], 2}
*)


It works fine if a different symbol is used:

c[x] /. First@sol /. x -> 0.5
(*  2.  *)
`