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) \\ \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)$.

 

[blue_raver22]

There are a few different approaches one could take to numerically solve a PDE with a delta function boundary condition. One option is to discretize the domain of the PDE and use a finite difference method, such as the one implemented in MATLAB's pdepe function. In this approach, the delta function boundary condition can be approximated by a very narrow and tall spike at the corresponding boundary point. This spike can then be represented by a large value at that point, while the rest of the boundary values remain unchanged.

Another approach is to use a finite element method, which is also commonly used for solving PDEs numerically. In this method, the domain is divided into smaller subdomains, and the PDE is approximated by a set of basis functions within each subdomain. The delta function boundary condition can be incorporated by appropriately choosing the basis functions at the boundary point where the delta function is located.

It is also possible to use a spectral method, where the solution is approximated by a combination of trigonometric or polynomial functions. In this approach, the delta function boundary condition can be incorporated by considering the Fourier or Chebyshev coefficients of the solution at the boundary point.

Regardless of the numerical method chosen, it is important to ensure that the boundary conditions are properly incorporated into the discretization scheme. In the case of a delta function boundary condition, special attention should be paid to accurately representing the sharp spike at the boundary point. Additionally, the numerical solution should be checked for stability and convergence to ensure that it accurately captures the behavior of the PDE with the delta function boundary condition.