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

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:

HallsofIvy

Science Advisor
Homework Helper
41,639
837
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.
 
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.
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
 
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].
 

Physics Forums Values

We Value Quality
• Topics based on mainstream science
• Proper English grammar and spelling
We Value Civility
• Positive and compassionate attitudes
• Patience while debating
We Value Productivity
• Disciplined to remain on-topic
• Recognition of own weaknesses
• Solo and co-op problem solving
Top