# How to correctly use DSolve when the force is an impulse (dirac delta) and initial conditions are not zero

DSolve (and NDSolve) return different and unexpected solution to differential equation when the input is an impulse. This is because DiracDelta[t] at $t=0$ remains DiracDelta.

This question is two parts question: If it is a user error to use DiracDelta[t] in this example, should then software still return a solution? And the second part asks how to use DSolve or NDSolve in order to obtain the correct solution to a differential equation when the input is an impulse $\delta\left(t\right)$ as is commonly understood and used in engineering and mathematics problems (the dirac delta function).

Initial conditions used are not all zero in this example.

## Problem description

An example which is standard ODE is used that models a single degree of freedom second order differential equation for an underdamped system. When using $f(t)=\delta\left( t\right)$ as the RHS of the differential equation, the answer returned is the same as the answer when $f\left( t\right) =0$. The user might go on and use the solution without knowing that $\delta\left( t\right)$ effect was not used since no warning is issued.

Here is an example solving $my^{\prime\prime}\left( t\right) +cy^{\prime }\left( t\right) +k=f\left( t\right)$. The solution was obtained for $f\left( t\right) =\delta\left( t\right)$ and $f\left( t\right) =0$ and for the same initial conditions $y\left( 0\right) =1,y^{\prime}\left( 0\right) =0$

The answers returned are the same as can be seen by the plots, and also the analytical answers are the same as well.

plot[title_, sol_] :=
Plot[sol, {t, 0, 10}, PlotRange -> All, Frame -> True,
FrameLabel -> {{y[t], None}, {Row[{t, " (sec)"}], title}},
GridLines -> Automatic, ImageSize -> 300];
eq1 = m y''[t] + c y'[t] + k y[t] == DiracDelta[t];
eq2 = m y''[t] + c y'[t] + k y[t] == 0;
parms = {m -> 1, c -> .7, k -> 5};
sol1 = First@DSolve[{eq1 /. parms, y == 1, y' == 0}, y[t], t];
sol2 = First@DSolve[{eq2 /. parms, y == 1, y' == 0}, y[t], t];
Grid[{
{plot["DiracDelta[t]", y[t] /. sol1],
plot["zero RHS", y[t] /. sol2]},
{Simplify[y[t] /. sol1], Simplify[y[t] /. sol2]}
}, Frame -> All] The analytical answer are the same. When impulse is input, the solution given is

E^(-0.35 t) (1. Cos[2.20851 t] - 0.294317 Sin[2.20851 t] +
0.452795 HeavisideTheta[t] Sin[2.20851 t])


and when the input is zero, the solution given is

E^(-0.35 t) (1. Cos[2.20851 t] + 0.158478 Sin[2.20851 t])


And they are the same since HeavisideTheta[t]=1 for $t>0$

The correct analytical solution for this differential equation derived by hand is different. And a plot of this solution is compared with the above solutions by Mathematica.

Taking the case of underdamped system, the analytical solution is $$u\left( t\right) =e^{-\xi\omega t}\left( u(0)\cos\omega_{d}t+\frac {u^{\prime}(0)+u(0)\xi\omega}{\omega_{d}}\sin\omega_{d}t\right) +e^{-\xi\omega t}\left( \frac{1}{m\omega_{d}}\sin\omega_{d}t\right)$$

where $\zeta=\frac{c}{2\omega m},\omega_{d}=\omega\sqrt{1-\zeta^{2}}% ,\omega=\sqrt{\frac{k}{m}}$, $u(0)$ and $u^{\prime}(0)$ are the initial conditions.

For the parameters in this example, the solution is

parms = {m -> 1, c -> .7, k -> 5};
w = Sqrt[k/m];
z = c/(2 m w);
wd = w Sqrt[1 - z^2];
analytical =
Exp[-z w t] (u0 Cos[wd t] + (v0 + (u0  z w))/wd Sin[wd t] +1/(m wd ) Sin[wd t]);
analytical /. parms /. {u0 -> 1, v0 -> 0}

(* E^(-0.35 t) (Cos[2.20851 t] + 0.611273 Sin[2.20851 t]) *)


Compare the above to the solution by DSolve which is given above, here it is again:

 (* E^(-0.35 t) (1. Cos[2.20851 t] + 0.158478 Sin[2.20851 t]) *)


Now plotting the above hand derived solution to compare:

Grid[{
{plot["DiracDelta[t]", y[t] /. sol1],
plot["drived", analytical /. parms /. {u0 -> 1, v0 -> 0}]}
}, Frame -> All] The difference can be seen clearly around $t=0$. A unit impulse has the effect of giving the system an initial velocity of $1/m$. In otherwords, a system with an impulse as input, can be replaced by an equivalent free input system by adding to its initial velocity an amount of $1/m$. This is why the initial velocity can not be zero in the solution. The slope of the displacement curve generated by the derived solution shows this, while the solution by DSolve used the zero initial velocity. This is all because DiracDelta[t] was basically not used.

ps. If I have to post the derivation of the solution done by hand, I can do this as well.

• Related? stackoverflow.com/q/14178369/353410 :) – Dr. belisarius Jan 8 '13 at 2:07
• Closely related:Differential equations with jump conditions – Jens Jan 8 '13 at 7:49
• The solution from DSolve is correct. Your initial conditions at $t=0$ prevent integration (from 0- to 0+) through the Dirac delta, i.e. there is no impulsive force seen by the system. However, if you move your initial conditions to $t<0$ then the Dirac delta does cause an impulsive force, which you can see as a discontinuity at $t=0$ when you plot $y'(t)$. – Stephen Luttrell Jan 8 '13 at 12:13

### What's going on

The real issue here is in how one should treat the Dirac delta in solving differential equations. You can check by integrating your equation over a small vicinity of zero that the presence of DiracDelta leads to a jump in derivative of your equation:

$$\int_{-eps}^{eps}(m y^{\prime\prime}(t)+c y^{\prime}(t)+k) dt == \int_{-eps}^{eps} \delta(t)dt$$

You can see that last 2 terms on the left vanish in the limit $$eps \rightarrow 0$$ - the term with k trivially, and the term with first derivative vanishes because the function itself is continuous at 0. However, the term with the second derivative does not vanish, but results in the jump of the first derivative:

$$m(y^{\prime}(+0)-y^{\prime}(-0)) = 1$$

where I used that the integral over the delta-function on the r.h.s.gives 1.

Physically this is correct, since the delta-function on the r.h.s. represents sudden infinite force acting during infinitely small time but giving the system non-zero change in momentum (velocity).

### How to get the correct solution in Mathematica

How does this look on the level of code: let us define a small epsilon and solve the equation for t<-eps and for t>eps separately:

eps=1.*10^-6

DSolve[{m y''[t]+c y'[t]+k y[t]==0/.parms,y[-eps]==1,y'[-eps]==0},y[t],t]

(* {{y[t]->1. E^(-0.35 t) (1. Cos[2.20851 t]+0.158476 Sin[2.20851 t])}} *)

DSolve[{m y''[t]+c y'[t]+k y[t]==0/.parms,y[eps]==1,y'[eps]==1},y[t],t]

(* {{y[t]->0.999999 E^(-0.35 t) (1. Cos[2.20851 t]+0.611276 Sin[2.20851 t])}} *)


where I used that m=1 and therefore the jump in the derivative is also 1.

While it would be nice if DSolve could handle such discontinuities automatically (perhaps it does and I just don't know how to specify this), my advice for now would be to resort to solving your equations separately in different domains, and glue the solutions together according to the boundary conditions induced by delta-function(s).

### Numerical approach

If all you need is a numerical solution, you can define an approximation for the delta-function, which would be centered at time eps and have a duration dt, as follows:

ClearAll[delta];
delta[t_, eps_, dt_] :=
1/dt*HeavisideTheta[t - eps]*HeavisideTheta[eps + dt - t]


You can now solve the equation without making any modifications to the initial conditions, as follows:

eps = 10^-6;
dt = 10^-7;
nsol =
NDSolve[{m y''[t] + c y'[t] + k y[t] == delta[t, eps, dt] /. parms,
y == 1, y' == 0}, y[t], {t, 0, 2 Pi}]


There is no ambiguity in this case, since the impulse occurs very shortly after zero time, but still not exactly at zero. You can check that the resulting solution is approximately correct. You can adjust the accuracy of this approximation by making eps and dt smaller.

### Conclusions

So, to reiterate: the answer produced by DSolve is not literally wrong, given the initial conditions you specified. It is the conditions which are ambiguous, since the correct solution has a first derivative that is discontinuous at zero (so that the first derivative can not be zero both at -0 and at +0). Solving in the two domains separately and explicitly taking into account this discontinuity leads to a correct solution. Of course, it would be nice if DSolve could handle such cases more autmoatically, but presumably it currently can not.

• @Nasser The problem is, in cases such as this one would need a special syntax. You tell NDSolve that the derivative must be zero at zero, but should it be zero before or after the impulse? This is ambiguous. More generally, delta-function is not really a normal function, it is a distribution. It has a clear meaning only together with the integration procedure - it is basically a kernel of an integral operator. But a step from a differential equation to integration is generally non-trivial, so I would be surprised if NDSolve managed to support delta-functions generally, even in principle. – Leonid Shifrin Jan 8 '13 at 3:31
• @Nasser In fact, usually differential equations involving delta-functions have a special status - usually those are equations for Green's functions of some differential operator, which are used then to solve the corresponding equations with some non-trivial r.h.s. So, I would not be so hard on DSolve regarding them, they are special. – Leonid Shifrin Jan 8 '13 at 3:36
• @Nasser I think this situation is an unfortunate consequence of popularization of delta-functions etc. which led to over-simplification of these concepts. Mathematica's documentation AFAIK did not state that DSolve can work with distributions (in the sense of general support), and I don't see why this needs to be mentioned explicitly. The real correct way to solve this problem is to consider an impulse of a finite time duration dt (pick any delta-function representation you like), which happens at time eps. When the solution is obtained, take dt to zero and eps to zero, in that ... – Leonid Shifrin Jan 8 '13 at 3:43
• @Nasser ... order. This should give you the correct solution, and in such a case you should be able to use NDSolve or DSolve without problems. But of course, this is a more time-consuming way to do it, so textbooks take short-cuts, and so do students. The bad thing here is that such short-cuts are detrimental to the understanding of things. In a way, I don't think it is bad that Mathematica does not fully automate this - delta-functions are subtle enough that one should always pay close personal attention to equations involving them - at least this is my view on it. – Leonid Shifrin Jan 8 '13 at 3:47

A straight forward way that MMa seems to do correctly is thus. Equation with DiracDelta with the impulse offset by a very small t0.

eq1 = m y''[t] + c y'[t] + k y[t] == DiracDelta[t - t0]


and the homogeneous eq.

eq2 = m y''[t] + c y'[t] + k y[t] == 0

parms = {m -> 1, c -> .7, k -> 5, t0 -> 10^-6}

DSolve[{eq1 /. parms, y == 1, y' == 0}, y[t], t] // Flatten;

y1[t_] = y[t] /. %;

DSolve[{eq2 /. parms, y == 1, y' == 0}, y[t], t] // Flatten;

y2[t_] = y[t] /. %;

Plot[{y1[t], y2[t]}, {t, 0, 10}, PlotRange -> All] To show that it honors the conditions at t = 0

Plot[{y1[t], y2[t]}, {t, 0, .000002}, PlotRange -> All] We can't have the DiracDelta spike at t == 0 and also specify y' == 0. This incompatibility is undoubtedly why DSolve did not do so well. I presume that since the DiracDelta fn in this case is part of the diffeq rather than a boundary or initial condition, there is no warning about inconsistent conditions. DSolve appears to do well in this case with consistent conditions, however.

It should be noted that t0 need not be near t = 0. Interesting solutions also occur if the impulse is at a later time.