I'm trying to solve the following equation (in Mathematica 10):

mysol = 
    {D[u[t, x], t] == 
       lambda*D[u[t, x], x, x] + p*(1 + Cos[2*omega*t])*DiracDelta[x],
     u[t, -b] == T0, u[0, x] == T0, Derivative[0, 1][u][t, a] == 0}, 
   {t, 0, 0.1}, {x, -b, a},
   Method -> {"PDEDiscretization" -> "FiniteElement"}]

with the constants:

T0 = 295; a = 100*10^(-9); b = 1*10^(-3); lambda = 0.5; p = 1*10^8; omega = 1000.

However the result I get is constant everywhere and I don't think it 'sees' the DiracDelta function. I've tried changing the method to the Finite Element method (as above) and replacing the DiracDelta function with D[HeavisideTheta[x],x] (as suggested elsewhere).

Any help would be greatly appreciated.

  • 4
    $\begingroup$ You can't use Dirac delta functions in numerical PDE solving routines directly, since Dirac delta function is not a normal function, but rather a distribution (kernel of an integral operator). What one usually does is to integrate the equation involving Dirac delta function in the small vicinity of its localization point, to get the boundary conditions for the solution to the left and to the right of it. Then, you solve separately on the left and on the right, and connect the solutions using those boundary conditions. Have a look at how Dirac-delta potential is solved in Quantum Mechanics. $\endgroup$ Oct 2, 2014 at 10:47
  • $\begingroup$ @LeonidShifrin Thanks- that's really helpful. Ultimately I need to extend this to 2D though, i.e. the Dirac Delta function becomes a line of points (with a finite width). Is there any way to deal with this? $\endgroup$
    – Susan
    Oct 2, 2014 at 12:14
  • 1
    $\begingroup$ It is not clear what you mean by finite width. Whenever you hit a delta-function, in any number of dimensions, it usually means that you have to divide your volume into regions, separated by the support of the delta-function, and then expect certain jumps in your function and / or derivatives (usually derivatives), so you'd need to solve separately and than match the integration constants, taking into account the contribution of the delta - function when you cross the boundary. This is as much as I can say on the general grounds. $\endgroup$ Oct 2, 2014 at 12:24
  • 2
    $\begingroup$ My final advice to you: study some examples of how Delta-functions are dealt with in the context of ODEs, first in 1D, then for PDEs in higher dimensions. Then, understand your full problem, and formulate it as a ODE (or PDE). Your question in its current stage has to do more with math than Mathematica per se. Make sure you get rid of Delta-function (translate its presence to certain jumps in boundary conditions) before feeding your equations to numerical solvers. $\endgroup$ Oct 2, 2014 at 12:45
  • 1
    $\begingroup$ Another possible way out here would be to approximate delta-function using some of the finite-width representations (many are well-known, all assume some limiting procedure of certain width parameter going to zero), then solve normally, then take a numerical limit of the width parameter going to zero. $\endgroup$ Oct 2, 2014 at 13:12


Browse other questions tagged or ask your own question.