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

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[0] == 1, y'[0] == 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.


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[0] == 1, y'[0] == 0}, y[t], t] // Flatten;

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

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

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

Plot[{y1[t], y2[t]}, {t, 0, 10}, PlotRange -> All]

enter image description here

To show that it honors the conditions at t = 0

Plot[{y1[t], y2[t]}, {t, 0, .000002}, PlotRange -> All]

enter image description here

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.