Differential equations with jump conditions

Suppose I want to solve an ODE with a DiracDelta source term. In the following example, DSolve does it correctly while NDSolve (understandably) misses the Delta and provides the solution corresponding to a zero right hand side:

ex = y /. First[DSolve[{y'[x] + y[x] == DiracDelta[x - 1], y == 1}, y, x]];
nu = y /. First[NDSolve[{y'[x] + y[x] == DiracDelta[x - 1] , y == 1}, y, {x, 0.1, 2}]];
Plot[ex[t], {t, 0.9, 1.1}]
Plot[nu[t], {t, 0.9, 1.1}]

Plot of DSolve solution Plot of NDSolve Solution Is there anyway around this using NDSolve (without converting the system to integral form)? I'd be happy to lend NDSolve a hand, say by telling it to look for a piecewise continuous solution ({x,0.1,1,2}) and provide jump conditions (in the above example y(1+eps) - y(1-eps) == 1).

(DSolve does not solve the real equations I am really interested in... they are multidimensional, nonlinear, and a bit too messy to replicate here).

• Which version are you using ? – b.gates.you.know.what Dec 1 '12 at 22:38
• version 8.(0.4.0) – Eric Dec 1 '12 at 22:45
• If possible, consider switching to Mathematica 9, there are capabilites to handle such problems in a seamless way. Take a look at this wolfram.com/mathematica/new-in-9/… – Artes Dec 1 '12 at 23:05

If you can't switch to Mathematica 9, then for numerical purposes it may be good enough to approximate the delta function by a strongly peaked, normalized function. The normal distribution is one that comes to mind:

nu = With[{ε = .0001},
y /. First[
NDSolve[{y'[x] + y[x] ==
PDF[NormalDistribution[1, ε], x], y == 1},
y, {x, 0.1, 2}, MaxStepSize -> ε,
MaxSteps -> Infinity]]];

Plot[nu[t], {t, 0.9, 1.1}] • +1. Was unaware of this ability of NDSolve, as well as your answer. – Leonid Shifrin Jan 8 '13 at 11:41