# Issues with modeling pulses in a very simple system of DAEs

Bug introduced in 10.4.1 or earlier and fixed in 11.0.0

Bug has been confirmed by WRI [Case:3594387]:

It does appear that the NDSolve function is not behaving properly in this case and an incident report has been created with the information you provided.

Some time ago I had asked this question about evaluation difficulties using Euler Integration to solve a system of ODE where discrete pulses occur.

While I have now abandoned Euler integration and can thus make use of DiracDelta to model pulses trying to introduce algebraic equations seems to pose a problem:

## Modeling an account with deposits

Consisder the most simple problem in Economics and Business, e.g. modeling a bank account with interest:

• $x(t)$ is to designate the amount of money in the bank account; at the beginning there will be 100 units of money in the account ($x(0) = 100$)
• there will be interest paid at a fractional rate ot $r = 0.05$ per unit of time
• There will be a deposit of +10 units at times $t_1 = 3$ and $t_2 = 6$
• The account is to be simulated for 10 periods

Setting this up with NDSolve is rather straight forward and works fine:

deposit = Function[t,
Total @ {
10. DiracDelta[ t - 3 ],
10. DiracDelta[ t - 6 ]
}
];

sim = First @ NDSolve[
{
x[0] == 100.,
x'[t] == 0.05 x[t] + deposit[t]
},
x,
{ t, 0, 10 }
];

Plot[ Evaluate @ ( x[t] /. sim ), { t, 0, 10 }]


## ... but will not work as DAEs (with a simple output function)

But what if this is set up as a differential algebraic equations model and there is to be some kind of system output as the engineers do in their control system framework (cf. NonlinearStateSpaceModel)?

More precisely: What happens if there is some function of the stock $y(t) = g(x(t))$ or even more simple some function of time $y(t) = g(t)$?

simDAE = First @ NDSolve[
{
x[0] == 100.,
x'[t] == 0.05 x[t] + deposit[t],

y[0] == 1.,
y[t] == 1. (* so our g() is simply a constant *)
},
{
x, y
},
{ t, 0, 10 }
];

Plot[ Evaluate @ ( x[t] /. simDAE ), {t, 0, 10} ]


Surprisingly this simple system is not simulated correctly begging the question:

What is going on here? Why are the pulses evaluated with the wrong sign?

UPDATE

Using another Method as has been suggested by MichaelE2 below will work for the very simple case so far, but unfortunately not for the more general case, e.g.

simDAE2 = First @ NDSolve[
{
(* modeling the system *)
x[0] == 100.,
x'[t] == 0.05 x[t] + deposit[t],

(* modeling the system's output *)
y[0] == x[0],
y[t] == x[t]
},
{ x, y},
{ t, 0, 10 },
Method -> { "EquationSimplification" -> "MassMatrix" }
];

Row @ Map[
Plot[ Evaluate @ (#[t] /. simDAE2 ),
{t, 0, 10},
ImageSize -> {GoldenRatio 200, 200},
PlotLegends -> Placed[ ToString @ #, Below]] &,
{x, y}
]


• My guess is that it' because your y[t] is not a differential equation, which may lead NDSolve to trying silly things... You'll note that if you define y[0]==1. and y'[t]==0 (which keeps it constant as you wish) that everything is ok – Quantum_Oli Apr 28 '16 at 11:56
• @Quantum_Oli Yes, I indeed noted, that ODEs do not present a problem, but you should be able to model $y(t) = g( x(t), u(t),t)$ as it is done for control system modeling in Mathematica (which internally uses NDSolveafaik?). Using output functions $y(t)$ will work fine for ODE without pulses btw. – gwr Apr 28 '16 at 12:09
• I updated the OP to account for the workaround proposed by @MichaelE2 which unfortunately does not seem robust. – gwr Apr 28 '16 at 12:55
• The issue has been confirmed as a bug by WRI as of today [Case:3594387]. – gwr Jun 8 '16 at 10:22
• While I have not gotten any other notice than that an incident report has been filed, the issue apparently has been fixed. At least the orignial example given now works fine in Version 11.0.0. – gwr Sep 5 '16 at 21:01

Below is a workaround for the simple case. The OP can say whether it works more general. I haven't quite tracked down yet why the system is set up incorrectly with the default Method -> {"EquationSimplification" -> "Solve"} and with Method -> {"EquationSimplification" -> "Residual"}. But it works in this case with Method -> {"EquationSimplification" -> "MassMatrix"}.

simDAE = First@
NDSolve[{x[0] == 100., x'[t] == 0.05 x[t] + deposit[t], y[0] == 1.,
y[t] == 1. }, {x, y}, {t, 0, 10},
Method -> {"EquationSimplification" -> "MassMatrix"}];

Plot[Evaluate@(x[t] /. simDAE), {t, 0, 10}]


Update: Response to comment

Again, I can only present a potential workaround at this point. Constructing a WhenEvent[] seems to work better than the automatic processing of DiracDelta[] in this case. For what it's worth, here's a function to convert deposit to a sequence of events, but it's probably easier to use its last line to construct the sequence directly from the times and amounts.

ClearAll[diracToEvent];
diracToEvent[depositFn_, x_, t_, scale_: 1/2] :=
Module[{
times =
Union @@ Cases[depositFn, DiracDelta[e_] :> (t /. Solve[e == 0, t]), Infinity],
dt,
amounts},
With[{xt = If[MatchQ[x, _[t]], x, x[t]]},
dt = Min@Differences@times;
amounts = Integrate[depositFn, {t, # - dt*scale, # + dt*scale}] & /@ times;
WhenEvent[t > #1, xt -> xt + #2] &,
{times, amounts}
]
]
];

simDAE = First@NDSolve[{
x[0] == 100., x'[t] == 0.05 x[t], y[t] == 25 Log[x[t]] ,
diracToEvent[deposit[t], x, t]},
{x, y}, {t, 0, 10}
];

Plot[Evaluate@({x[t], y[t]} /. simDAE), {t, 0, 10}]


• That works for the simple case I gave, but if you try y[0] == x[0], y[t] == x[t] it does not look very reassuring. In a more complex model I am also noting that "EquationSimplification" -> "MassMatrix" unfortunately does not reliably fix the problems with DiracDelta pulses. – gwr Apr 28 '16 at 12:36
• A more robust way to convert an expression with linear DiracDelta[] terms might be amounts = Limit[SeriesCoefficient[depositFn, {DiracDelta[t - #], 0, 1}], t -> #] & /@ times. My original thought was that it would be fun to try to make a general converter; but there's probably no interest, and the fun lasts only so long. :) – Michael E2 Apr 28 '16 at 14:11
• Thank you, the WhenEvent solution works nicely. Would you consider the described behavior for DiracDelta a bug? I may underestimate the difficulties in DAEs, but it seems simple enough here so that DAE vs. ODE should not make a difference maybe? – gwr May 2 '16 at 10:56
• @gwr I think it's worth reporting to WRI. I don't know enough (not anything, really) about the internal code to say for sure whether this is a bug or a limitation -- that is to say, whether it's a programming mistake or an edge-case of a general algorithm which turns out to be difficult to handle. The problem seems a common enough type that they might want to fix it. – Michael E2 May 2 '16 at 12:07
• I had reported the issue [Case:3594387] and it has today been confirmed as an issue within the NDSolve function. An incident report has been created so there is hope for the future which should not die. :) – gwr Jun 8 '16 at 10:18