Wrong derivative in delay differential equation initial data returned by NDSolve

Question

I have an DDE:

sol = NDSolve[{x'[t] + x[t - 1] == 0, x[t /; t <= 0] == t^2}, x, {t, -1, 5}]

And result (take the curve on the time interval $-1<t<0$)

Plot[Evaluate[{x[t], x'[t]} /. sol], {t, -1, 0}, PlotRange -> All]

The following is not clear to me:

We have $x(t)=t^2$ and why $x'(t)=-2$, but $x'(t)≠2t$ ?

The history is stored in an Experimental`NumericalFunction. The function is defined to be t^2, but there seems to be a bug with its derivative. See this: ifn = x /. First@sol; nf = ifn[[4, 1, 2]]; nf["Jacobian"[-3]]. Whatever number is used instead of -3, you get -2. for the derivative ("Jacobian"). Please report it to WRI. — Michael E2, Jul 12 '22 at 04:22
@MichaelE2 Really curious. So this is really a bug? At first I thought that this is a feature that the delay introduces. — dtn, Jul 12 '22 at 04:23
The numerical function should work as a regular function, don't you think? It shouldn't have a random value for the derivative. — Michael E2, Jul 12 '22 at 04:26
@dtn If we define solution at t<=0 as x[t /; t <= 0] == t^2 then we can integrate delay equation at t>0 only. In your example you try to integrate at {t, -1, 5}. — Alex Trounev, Jul 12 '22 at 05:05
@AlexTrounev I reason based on some observations. Numerical analysis shows that on the time interval $-1<t<0$ $x(t)$ behaves in accordance with the given history (in our case it is $t^2$). At the zero moment of time, the time domain familiar to us for numerical calculations begins, and one of the accepted numerical algorithms begins to work. By the zero moment of time, both the initial value of the calculated function and its delayed value are available to us. — dtn, Jul 12 '22 at 06:31
@AlexTrounev If the numerical method is correctly applied, then the condition laid down in the DDE will be fulfilled. This refers to $t>0$. As for the time $-1<t<0$, some kind of contradiction is intuitively felt there, that it is impossible to fulfill several conditions $x'(t)=-x(t-1)$ and $x'(t)=x(t)$ at the same time. Those, when $x'(t)=-x(t-1)$, then $x'(t)≠x(t)$, and so on. All conditions were satisfied with exponential functions, probably. Therefore, with $-1<t<0$, an assumption is made about the discrepancy between the trajectory of the desired function and its derivative. — dtn, Jul 12 '22 at 06:34
@AlexTrounev I reproduced the result from Mathematica in Excel, based on these considerations (with Euler Finite Difference Method): https://ibb.co/2Zkfgjz — dtn, Jul 12 '22 at 06:37
@dtn Please, read attentively theorem about uniqueness solution for DDE, for example, here https://www.researchgate.net/publication/309501592_Existence_and_uniqueness_of_the_solution_of_delay_differential_equations — Alex Trounev, Jul 12 '22 at 10:40
The same issue with sol = NDSolve[{x'[t] + x[t - 1] == 0, x[t /; t <= 0] == Sin[t]}, x, {t, -1, 5}]. — user64494, Jul 12 '22 at 11:53
@AlexTrounev You might be interested that the docs show an example like the OP's (but with x''[t] instead of x'[t]): sol = NDSolve[{x''[t] + x[t - 1] == 0, x[t /; t <= 0] == t^2}, x, {t, -1, 5}] (first example here). It, too, produces a numerical function with a buggy nf["Jacobian"[-1]]. I don't think specifying the domain to include the history means NDSolve will integrate the DE. It simply incorporates the initial history and integrates only from t == 0 forward. — Michael E2, Jul 12 '22 at 15:10
@MichaelE2 Yes, they have example with x[t /; t <= 0] == t^2 and {t, -1, 5} solved with NDSolve, but that is incorrect. In DDE theory we can't solve equation at t<=0 since solution is defined at t<=0 as an initial condition. — Alex Trounev, Jul 12 '22 at 16:32
@AlexTrounev I understand. NDSolve does not solve the DDE for t <= 0 even if you specify {t, -1, 5}. While it is incorrect to solve the DDE for t <= 0, NDSolve is not doing the incorrect thing. So there is nothing incorrect with specifying {t, -1, 5}. It simply uses the initial history for t <= 0, which is correct, yes? — Michael E2, Jul 12 '22 at 16:45

Michael E2 · Accepted Answer · 2022-07-12T20:11:27.903

I figured out part of the problem, and it no longer seems random. My first answer-comment was premature, but I was confused and thought I had gotten as far as I could. The problem occurs as described below. I don't know why it messes up the numerical function's "Jacobian".

Inside the InterpolatingFunction solution, there is an Experimental`NumericalFunction that encodes the history of the DDE; let's call it nf[]. The problem occurs when nf["Jacobian"[...]] is evaluated with whatever numeric argument for ... while t has been given a value.

Once this happens, then the Jacobian's value for x'[...] is fixed. You can rerun NDSolve to create a new solution with a new numerical function.

On the other hand, if the Jacobian is first evaluated when t has no value, then it will evaluate correctly from that point on, even if t is given a value afterwards.

Note that the numerical function is used to evaluate the derivative x'[t] only for values of t in the initial history, that is, for t < 0 in the OP's problem. For higher order DDEs, the Jacobian is not evaluated for derivatives of order less than the order of the DDE. This problem arises in higher order DDEs only when evaluating derivatives of the order of the DDE or higher.

Evidence for the above:

We start over with a new NDSolve every time to get a new numerical function in the solution.

Clear[x, t];
sol = NDSolve[{x'[t] + x[t - 1] == 0, x[t /; t <= 0] == t^2}, 
   x, {t, -1, 5}];
t = -0.44;
x'[-0.6] /. First@sol          (* t has a value *)
x'[{-0.5, -0.33}] /. First@sol (* all derivatives = x'[-0.44] *)
t =.
(*
  -0.88
  {-0.88, -0.88}
*)
Clear[x, t];
sol = NDSolve[{x'[t] + x[t - 1] == 0, x[t /; t <= 0] == t^2}, 
   x, {t, -1, 5}];
x'[-0.6] /. First@sol; (* t has no value )
t = -0.44;
x'[-0.6] /. First@sol  ( the value of t is irrelevant )
x'[{-0.5, -0.33}] /. First@sol ( & derivatives are correct )
t =.
(
  -1.2
  {-1., -0.66}
*)
Clear[x, t];
sol = NDSolve[{x'[t] + x[t - 1] == 0, x[t /; t <= 0] == t^2}, 
   x, {t, -1, 5}];
Block[{t = 0.4},        (* positive t but... )
 x'[-0.6] /. First@sol] ( evaluate x' at neg. number )
x'[{-0.5, -0.33}] /. First@sol
(
  0.8
  {0.8, 0.8}
*)

Another test: Using a different variable than t in Plot avoids the issue (unless the x' has been evaluated at a negative number when t had a value already, of course):

Clear[x, t];
sol = NDSolve[{x'[t] + x[t - 1] == 0, x[t /; t <= 0] == t^2}, 
   x, {t, -1, 5}];
Plot[Evaluate[{x[tt], x'[tt]} /. First@sol], {tt, -1, 2}, 
 PlotRange -> All]

I'm deleting most of the original answer, because the above is more accurate and clearer. However, there is a workaround I had proposed that I cannot fully explain in terms of what was said above, but I have a guess. The following produces a solution that does not have the bug:

sol2 = NDSolve[
 {x'[t] + $MinMachineNumber*x'[t - 1] + x[t - 1] == 0, 
  x[t /; t <= 0] == t^2}, x, {t, -1, 5}]

If the first thing you do is plot it in terms of t, you get the same graph as above (not shown below):

Plot[Evaluate[{x[t], x'[t]} /. First@sol2], {t, -1, 2}, 
 PlotRange -> All]

If we run one of the tests above, we see we get different values from the two numerical functions even though according to SameQ they are identical. Something must happen internally, which I cannot discover, that makes the difference. Perhaps the Jacobian of the numerical function is evaluated to determine x[t - 1] during the integration of the DDE. It seems plausible and an internal Trace shows that a numerical function like the one in the solution occurs (in NDSolve`DDEFunction at each integration step). However, if it's evaluated, Trace does not reveal it. (Trace cannot show all steps.)

Clear[x, t];
sol = NDSolve[{x'[t] + x[t - 1] == 0, x[t /; t <= 0] == t^2}, 
   x, {t, -1, 5}];
Block[{t = -0.4}, x'[t] /. First@sol];
x'[{-0.5, -0.33}] /. First@sol
ifn = x /. First@sol;
nf = ifn[[4, 1, 2]];
nf["Jacobian"[#]] & /@ {-0.5, -0.33}
(*
  {-0.8, -0.8}
  {{-0.8}, {-0.8}}
*)
Clear[x, t];
sol2 = NDSolve[{x'[t] + $MinMachineNumberx'[t - 1] + x[t - 1] == 0, 
    x[t /; t <= 0] == t^2}, x, {t, -1, 5}];
Block[{t = -0.4}, x'[t] /. First@sol2];
x'[{-0.5, -0.33}] /. First@sol2
ifn2 = x /. First@sol2;
nf2 = ifn2[[4, 1, 2]];
nf2["Jacobian"[#]] & /@ {-0.5, -0.33}
(
  {-1., -0.66}
  {{-1.}, {-0.66}}
*)
nf === nf2
(*
  True
*)

Do you mean sol1 = NDSolve[{x'[t] + 0. x'[t - 1] + x[t - 1] == 0, x[t /; t <= 0] == t^2}, x, {t, -1, 5}]; or sol1 = NDSolve[{x''[t] + 0. x'[t - 1] + x[t - 1] == 0, x[t /; t <= 0] == t^2}, x, {t, -1, 5}];? The send order or the first order? — user64494, Jul 12 '22 at 16:03
For me, in version 13.1 on Windows 10 sol1 = NDSolve[{x'[t] +0.*x'[t]++ x[t - 1] == 0, x[t /; t <= 0] == t^2}, x, {t, -1, 5}];Plot[Evaluate[{x[t], x'[t]} /. sol1], {t, -1, 2}, PlotRange -> All] produces the same phenomenon as in the question. The executed nb/ its pdf file on demand. — user64494, Jul 12 '22 at 16:21
Some reports: Take a look. I do not vouch for the absolute correctness, and perhaps there are flaws in my reasoning, but nevertheless. Based on the "historical" properties of DDE, we can assume that on the time interval (taking into account the selected delay, in our case $-1$) the history is $x(t)=t^2$. This means that on the time interval $-1≤t≤0$ $x'(t)=\frac{d}{dt}t^2$. — dtn, Jul 12 '22 at 17:33
Let's try to match the trajectory of a variable and its derivative using a piecewise function. As for the time interval $t>0$, this is the area of customary numerical analysis. I.e. $x'(t)+x(t-1)=0$ must be executed on the time interval $t>0$. Let's pay attention to the initial value of the $x(t)$ variable. It must equal the value of her "history" at time $t=-1$. Thus: solp = NDSolve[{x'[t] == Piecewise[{{D[t^2, t], -1 <= t <= 0}, {-x[t - 1], t > 0}}], x[-1] == (-1)^2}, x, {t, -1, 5}] And the result: Plot[Evaluate[{x[t], x'[t]} /. solp], {t, -1, 5}, PlotRange -> All, ImageSize -> Small] — dtn, Jul 12 '22 at 17:33
As I said, there may be flaws in the reasoning, but this approach can become the basis for a complete solution to the problem. Please, compare my solution with yours. — dtn, Jul 12 '22 at 17:35
@dtn The plots overlap exactly. It seems right. -- I guess it's wrong at one point only, x'[-1] /. solp. — Michael E2, Jul 12 '22 at 17:50
@dtn This fixes it: NDSolve[{x'[t] == Piecewise[{{D[t^2, t], -1 <= t <= 0}, {-x[t - 1], t > 0}}], x[t /; t <= -1] == t^2, x'[t /; t <= -1] == 2 t}, x, {t, -1, 5}] -- The same trick does not fix your original DDE, though. — Michael E2, Jul 12 '22 at 17:56
Indeed,sol = NDSolve[{x'[t] + $MinMachineNumber*x'[t - 1] + x[t - 1] == 0, x[t /; t <= 0] == t^2}, x, {t, -1, 5}]; is a tricky workaround. The edit in the answer $MinMachineNumber*x'[t - 1] instead of 0.*x'[t - 1] should be described. — user64494, Jul 12 '22 at 18:33
@dtn I figured out somewhat how the problem arises, though not fully why, and that leads to a better workaround, I think. Just evaluate the derivative at a negative number immediately after NDSolve: x'[-0.2] /. sol. Then the problem seems fixed in any further evaluations. — Michael E2, Jul 12 '22 at 20:15
I correctly understood that there is something inside the program that causes such glitches? Let's call it a small imperfection of the algorithm. And there are several ways to get around this problem, which we have discussed here? — dtn, Jul 13 '22 at 05:13

Wrong derivative in delay differential equation initial data returned by NDSolve

1 Answers1