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[0].

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.

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.

$\begingroup$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)$.$\endgroup$
– Stephen LuttrellJan 8 '13 at 12:13

2 Answers
2

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:

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:

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:

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.

$\begingroup$@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.$\endgroup$
– Leonid ShifrinJan 8 '13 at 3:31

$\begingroup$@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.$\endgroup$
– Leonid ShifrinJan 8 '13 at 3:36

$\begingroup$@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 ...$\endgroup$
– Leonid ShifrinJan 8 '13 at 3:43

$\begingroup$@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.$\endgroup$
– Leonid ShifrinJan 8 '13 at 3:47

We can't have the DiracDelta spike at t == 0 and also specify y'[0] == 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.

Mathematica is a registered trademark of Wolfram Research, Inc. While the mark is used herein with the limited permission of Wolfram Research, Stack Exchange and this site disclaim all affiliation therewith.