# Unsure solution of ODE with DiracDelta function

Could the following solution be treated as reliable?

DSolve[{D[x[t], t] == x[t]^3*DiracDelta[t - 2], x[1] == 1}, x[t], t]


DSolve::bvnul: For some branches of the general solution, the given boundary conditions lead to an empty solution.

{{x[t] -> 1/Sqrt[1 - 2 HeavisideTheta[-2 + t]]}}

I have in mind that the function x[t] takes complex values for $t>2.$

The general solution is produced with no warning:

DSolve[D[x[t], t] == x[t]^3*DiracDelta[t - 2], x[t], t]


{{x[t] -> -(1/(Sqrt[2] Sqrt[-C1 - HeavisideTheta[-2 + t]]))}, {x[t] -> 1/(Sqrt[2] Sqrt[-C1 - HeavisideTheta[-2 + t]])}}

Addition. Its numeric solution is quite different:

s = NDSolve[{D[x[t], t] == x[t]^3*DiracDelta[t - 2], x[1] == 1},x[t], {t, 1, 3}];
Plot[Evaluate[x[t] /. s], {t, 1, 3}, PlotRange -> All]


Which one of these solutions is more true?

• In the second, there are two $\pm$ branches; only one takes real values for some t when a condition is x[1] == 1l that leads to C[1] == -1/2. The second solution is empty for this initial condition, and this is what the first message is about. Next, not every differential equation has a solution that is valid on $(-\infty,+\infty)$, and this is just the case here. – corey979 Nov 30 '16 at 18:09
• @corey979: Thank you for your interest and feedback. The open question is: where is the obtained solution valid? – user64494 Nov 30 '16 at 18:12
• FunctionDomain[ 1/Sqrt[1 - 2 HeavisideTheta[-2 + t]] /. HeavisideTheta -> UnitStep, t] gives t < 2. Mathematica doesn't want to work with HeavisideTheta in FunctionDomain, hence the replacement to an equivalent UnitStep. – corey979 Nov 30 '16 at 18:22
• @corey979: I don't see a problem here: HeavisideTheta[1] outputs 1. – user64494 Nov 30 '16 at 18:37
• Ad Addition. Every piecewise (with regard to $t\lessgtr 2$) constant function is a solution to the ODE: $\delta(t-2)=0$ for $t\neq 2$, and $x'(t)$ of an arbitrary constant is equal to zero. $\delta(t)$ is not a function - it's a distribution, and common sense (i.e., treating it as a function) usually fails. The implementation of DiracDelta is connected with HeavisideTheta via HeavisideTheta'[t] == DiracDelta[t]. The internal algorithms of DSolve seem to rely on this relation. – corey979 Nov 30 '16 at 19:28

Converting DiracDelta[t-2] to Piecewise[{{DiracDelta[0], t == 2}}] helps.

DSolve[{D[x[t], t] == x[t]^3*Piecewise[{{DiracDelta[0], t == 2}}], x[1] == 1}, x[t], t]

(* {{x[t] -> 1}} *)


Okay so let me begin.

Mathematica itself does a semi-good job here. It find 2 solutions where one is not applicable for your startingcondition. So it finds the only solution which is

1/Sqrt[1 - 2 HeavisideTheta[-2 + t]]


But this only seems to work for $t<2$ and this is indeed right (we'll see later).

If we apply the Fourier- or Laplactransform we get what we deserve. (nothing we can work with because the $\delta$-function forces us to know $x(2)$)

$$x(t)=\theta(t-2)\cdot x(2)+C$$ Since 2 is a exact discontinuity we cannot use this. But we can use the good old approach of seperation of variables so we solve it by hand: $$\frac{\text{d}x(t)}{\text{d}t}=x^3(t)\cdot\delta(t-2)$$ $$\int\frac{1}{x^3(t)}\text{d}x(t)=\int\delta(t-2)\text{d}t$$ $$-\frac{1}{2x^2(t)}=\theta(t-2)+C$$ Now we would like to use $(...)^{-1}$ but this will limit us to the domain of $t>2$ but this is okay, Mathematica solved it already for $t<2$. We get: $$x(t)=\pm\sqrt{C-\frac{\theta(t-2)}{2}}$$ with your condition: $$x(t)=\pm\sqrt{1-\frac{\theta(t-2)}{2}}=const((2<t)\wedge(2>t))$$

So we look into the visualization:

Plot[{1/Sqrt[1-2 HeavisideTheta[-2+t]],Sqrt[1-HeavisideTheta[t-2]/2]},{t,-3,3},PlotRange->{{-3,3},{-5,5}},PlotStyle->{{Red},{Blue,Dashed}},PlotLegends->"Expressions"]


Oh look we have won! Our derived expression suprisingly also describes $t<2$. Let's confirm that:

equ=x'[t]==x[t]^3*DiracDelta[t-2];
form=Sqrt[1-HeavisideTheta[t-2]/2];
FullSimplify[equ/.{x[t]->form,x'[t]->D[form,t]},{t\[Element]Reals,t<2}]
FullSimplify[equ/.{x[t]->form,x'[t]->D[form,t]},{t\[Element]Reals,t>2}]
FullSimplify[equ/.{x[t]->form,x'[t]->D[form,t]},{t\[Element]Reals,t==2}]


True

True

-(DiracDelta[0]/(2 Sqrt[4-2 HeavisideTheta[0]]))==DiracDelta[0] (1-HeavisideTheta[0]/2)^(3/2)

Year. Exactly what we expected. Our formula desribes it perfectly EXCEPT in $t=2$.

• Note also my comment, and that $x(t)+c=\int x'(t)dt=\int x^3(t)\delta(t-2)dt=x^3(2)$, hence every piecewise constant function is a solution, e.g. $x(t)=-googol$ for $t<2$ ,and $x(t)=Knuth\, number$ for $t>2$. With the initial condition $x(1)=1$ the part for $t<2$ is $x(t)=1$, but for $t>2$ it still can be any number. – corey979 Nov 30 '16 at 20:15
• @JUlien Kluge: Thank you. As it was noticed by corey979, every piecewise (with regard to t≶2,t≶2) constant function is a solution to the ODE. – user64494 Nov 30 '16 at 20:15

Both the Mathematica answers in my question are not correct and should be treated as bugs. Here are arguments. Up to Rudin W. Functional Analysis, McGraw-Hill Book Company, NY,…,Toronto, 1973, Part 2, Ch. 6, Par. 6.15, the multiplication $f(t)\delta(t)$ is defined if the function $f(t)$ is infinitely differentiable. This implies the solution $x(t)$ must be continuous at $t=2$. Therefore, the only true solution of the ODE with the initial condition under consideration is a constant function $x(t)=1$. This is not only my personal opinion: I consulted my colleagues about this topic. In particular, one may contact Dr. Yu. Golovatyj at yu_holovaty@franko.lviv.ua (with his kind permission) for more details concerning this complex matter under consideration. Wolfram staff also obtained my report [CASE:3789146].

Mathematica's first (analytic) solution is dead on correct. It satisfies the differential equation and the boundary condition. When that is the case, you have the solution.

ode = D[x[t], t] == x[t]^3*DiracDelta[t - 2];

DSolve[{ode, x[1] == 1}, x[t], t] // Flatten
(*{x[t] -> 1/Sqrt[1 - 2*HeavisideTheta[t - 2]]}*)

x[t_] = x[t] /. %;


MMa thinks it satisfies the ode.

ode
(*True*)


and the bc

x[1]
(*1*)


Take the derivative by hand as a check. write x as

(1 - 2 HeavisideTheta[t - 2])^(-1/2)


and the derivative

(-1/2) (1 - 2 HeavisideTheta[t - 2])^(-3/2) (-2 DiracDelta[t - 2])
(*DiracDelta[t - 2]/(1 - 2*HeavisideTheta[t - 2])^(3/2)*)


We used the fact the derivative of the HeavisideTheta is the DiracDelta. Look at

x[t]^3 DiracDelta[t - 2]
(*DiracDelta[t - 2]/(1 - 2 HeavisideTheta[t - 2])^(3/2)*)

%%==%
(*True*)


Plot both real and imaginary parts of the solution

Looking at the ode, in this case the slope is 0 everywhere except at x = 2 and the solution matches that, but it does jump from pure real for x < 2 to pure imaginary for x > 2.
And no, x will not be continuous across the delta function. f[t] being continuous does not mean f[t] δ[t] is continuous.

The second (numerical) solution appears to be incorrect for values > 2. Even though numerical, the solution is evidently:

xnum[t_] = 1 + HeavisideTheta[t - 2]


check

ode /. x -> xnum

(*DiracDelta[t - 2] == DiracDelta[t - 2] (HeavisideTheta[t - 2] + 1)^3*)


Which in general is not true for all t.