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< -\epsilon$, = 1 if $-\epsilon< x< \epsilon$, = 0 if $x> \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< -\epsilon$, = 1 if $-\epsilon< x< \epsilon$, = 0 if $x> \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 < x < \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)$.