2
$\begingroup$

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

$\endgroup$
  • $\begingroup$ Which version are you using ? $\endgroup$ – b.gates.you.know.what Dec 1 '12 at 22:38
  • $\begingroup$ version 8.(0.4.0) $\endgroup$ – Eric Dec 1 '12 at 22:45
  • $\begingroup$ 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/… $\endgroup$ – Artes Dec 1 '12 at 23:05
5
$\begingroup$

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

$\endgroup$
  • $\begingroup$ +1. Was unaware of this ability of NDSolve, as well as your answer. $\endgroup$ – Leonid Shifrin Jan 8 '13 at 11:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.