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

[Only a Mirage]

I have a PDE of the following form:

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) \\<br /> \lim_{t\to s^+} f(t,x,y) = \delta (x-y)

Here k and c are real numbers and g, h 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?

 

[HallsofIvy]

You can approximate the delta function by a step function: \delta(x)= 0 if x&lt; -\epsilon, = 1 if -\epsilon&lt; x&lt; \epsilon, = 0 if x&gt; \epsilon for small \epsilon. How small \epsilon must be depends upon how small you are taking the step size in your numerical integration.

 

[Only a Mirage]

You can approximate the delta function by a step function: \delta(x)= 0 if x&lt; -\epsilon, = 1 if -\epsilon&lt; x&lt; \epsilon, = 0 if x&gt; \epsilon for small \epsilon. How small \epsilon must be depends upon how small you are taking the step size in your numerical integration.

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

 

[the_wolfman]

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 \int \delta(x) dx =1. So you should use
f(x) =\frac{ 1}{2\epsilon} for -\epsilon &lt; x &lt; \epsilon

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 f_i =\int w(x)\delta(x) \alpha_i(x) dx where \alpha_i(x) is your i-th basis function, w(x) is a weighting function, and f(x) =\sum_i f_i \alpha_i(x).