2

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[2] == 1}, y, x]];
nu = y /. First[NDSolve[{y'[x] + y[x] == DiracDelta[x - 1] , y[2] == 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 DSolve solution

Plot of NDSolve Solution

enter image description here

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

halirutan
  • 112,764
  • 7
  • 263
  • 474
Eric
  • 333
  • 1
  • 10

1 Answers1

6

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[2] == 1}, 
      y, {x, 0.1, 2}, MaxStepSize -> ε, 
      MaxSteps -> Infinity]]];

Plot[nu[t], {t, 0.9, 1.1}]

nu

Jens
  • 97,245
  • 7
  • 213
  • 499