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]

enter image description here

Which one of these solutions is more true?

  • $\begingroup$ 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. $\endgroup$
    – corey979
    Commented Nov 30, 2016 at 18:09
  • $\begingroup$ @corey979: Thank you for your interest and feedback. The open question is: where is the obtained solution valid? $\endgroup$
    – user64494
    Commented Nov 30, 2016 at 18:12
  • $\begingroup$ 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. $\endgroup$
    – corey979
    Commented Nov 30, 2016 at 18:22
  • $\begingroup$ @corey979: I don't see a problem here: HeavisideTheta[1] outputs 1. $\endgroup$
    – user64494
    Commented Nov 30, 2016 at 18:37
  • 1
    $\begingroup$ 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. $\endgroup$
    – corey979
    Commented Nov 30, 2016 at 19:28

4 Answers 4


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:




-(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$.

  • 1
    $\begingroup$ 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. $\endgroup$
    – corey979
    Commented Nov 30, 2016 at 20:15
  • $\begingroup$ @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. $\endgroup$
    – user64494
    Commented Nov 30, 2016 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 [email protected] (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.


and the bc


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)*)


Plot both real and imaginary parts of the solution

enter image description here

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]


ode /. x -> xnum

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

Which in general is not true for all t.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.