2

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

mysol = 
  NDSolve[
    {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}, 
   u,
   {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.

m_goldberg
  • 107,779
  • 16
  • 103
  • 257
Susan
  • 93
  • 7
  • 4
    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. – Leonid Shifrin Oct 02 '14 at 10:47
  • @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? – Susan Oct 02 '14 at 12:14
  • 1
    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. – Leonid Shifrin Oct 02 '14 at 12:24
  • @LeonidShifrin Sorry, I mean that I have two rectangular regions, with a line (the heater) sandwiched between them. However while the regions can be thought of as infinite in the y direction, the heater has a finite length in y (so there are some places where the two regions are in direct contact with each other, rather than being separated by the heater). Therefore in 1D the heater separates the 2 regions completely, but not in 2D. – Susan Oct 02 '14 at 12:39
  • 2
    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. – Leonid Shifrin Oct 02 '14 at 12:45
  • 1
    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. – Leonid Shifrin Oct 02 '14 at 13:12
  • @Susan, could you provide a geometry of your region in 2D? – user21 Oct 02 '14 at 15:58
  • Are you looking for DiscreteDelta? – bobthechemist Oct 02 '14 at 18:53

0 Answers0