0
$\begingroup$

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: enter image description here

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. >>
$\endgroup$
  • 1
    $\begingroup$ It appears the green line reaches 0 after t=1, causing a stiff system. $\endgroup$ – Feyre Oct 17 '16 at 10:00
2
$\begingroup$

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

Re NDSolve::ivcon:

Your initial conditions are inconsistent:

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[1][s][0] == 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[1][s][0], 0}}  *)

Re NDSolve::ndsz:

Mathematica graphics

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

Mathematica graphics

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 NDSolve`InitialHistory 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, 
   Experimental`NumericalFunction[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.  *)
$\endgroup$

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.