Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

How to numerically solve a PDE with delta function boundary condition?

  1. Jun 11, 2014 #1
    I have a PDE of the following form:

    [tex] f_t(t,x,y) = k f + g(x,y) f_x(t,x,y) + h(x,y) f_y(t,x,y) + c f_{yy}(t,x,y) \\
    \lim_{t\to s^+} f(t,x,y) = \delta (x-y)[/tex]

    Here [itex]k[/itex] and [itex]c [/itex] are real numbers and [itex]g, h[/itex] are (infinitely) smooth real-valued functions. I have been trying to learn how to do this numerically (was going to use MATLAB's pdepe function which uses a finite difference method, I believe), but I have no clue what to do regarding the boundary condition.

    How does one numerically integrate a PDE with dirac delta boundary condition?
    Last edited: Jun 11, 2014
  2. jcsd
  3. Jun 11, 2014 #2


    User Avatar
    Science Advisor

    You can approximate the delta function by a step function: [itex]\delta(x)= 0[/itex] if [itex]x< -\epsilon[/itex], = 1 if [itex]-\epsilon< x< \epsilon[/itex], = 0 if [itex]x> \epsilon[/itex] for small [itex]\epsilon[/itex]. How small [itex]\epsilon[/itex] must be depends upon how small you are taking the step size in your numerical integration.
  4. Jun 11, 2014 #3
    Oh okay, that makes sense. So that seems to solve that problem. Now I just need to find some software to do this in since apparently with the Matlab toolboxes I have I can only solve 1-dimensional PDEs
  5. Jun 13, 2014 #4
    There are a couple of different ways to represent a delta function initial condition numerically, but the representation should be consistent with the numerical method that you are using to solve the PDE.

    For instance using a step function as Hallsofivy suggested for finite difference methods. However there is an error.
    You want to preserve the integral [itex]\int \delta(x) dx =1[/itex]. So you should use
    [itex]f(x) =\frac{ 1}{2\epsilon}[/itex] for [itex]-\epsilon < x < \epsilon[/itex]

    If you're using spectral methods or finite element methods you should project the delta function onto each of your basis functions. The projection will look like [itex]f_i =\int w(x)\delta(x) \alpha_i(x) dx[/itex] where [itex] \alpha_i(x)[/itex] is your i-th basis function, [itex] w(x)[/itex] is a weighting function, and [itex]f(x) =\sum_i f_i \alpha_i(x)[/itex].
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook