Represent Dirac Delta function in Finite Difference method

Question

I recently solving $-\Delta u=\delta$ where $\delta$ is dirac delta function using FDM on 2 dimensional space.

Since dirac delta function is undefined at origin, and 0 elsewhere, I will use

$\delta(x)=\begin{cases}\delta_h\quad\text{at }(0,0)\\0\quad\text{elsewhere}\end{cases}$.

However, I don't know which value should I take for $\delta_h$. Options that I thinking are $\frac{3}{4\Delta x\Delta y}$(consider function shape as pyramid) and $\frac{1}{\Delta x\Delta y}$(consider function shape as cuboid)

Answer 1

First of all, for this problem polar coordinates are better suited than cartesian coordinates due to the radial symmetry of $\delta$. This is also argued in this notes which derive the fundamental solution

$$ u (x, y) = - \frac{1}{2\pi} \ln \bigg( \sqrt{x^2 + y^2} \bigg). $$

To simulate your problem, you need to impose some boundary conditions (which would find their way into the fundamental solution). The version of the fundamental solution stated above corresponds to homogeneous Dirichelt boundary condition on the unit circle: $u\bigg(\sqrt{x^2 + y^2} = 1\bigg) = 0$.

Coming back to your original question: I digged out some old cartesian FDM Matlab code of mine and ran it on the domain $[-1, 1]^2$ with Dirichlet boundary condition $0 = u(\pm 1, y) = u(x, \pm 1)$ for your different proposals. For the standard 5-point stencil $\delta_h = 1 / (\Delta x)^2$ worked (although the radial symmetry is clearly lost for the choice of my boundary conditions), see the code & picture below. It would be interesting to check also the third order FDM. Since this is also symmetric for the interior points and scaled with $1/(\Delta x)^2$, I suspect that again $\delta_h = 1 / (\Delta x)^2$ would work.

clear;
clc;
close all;

n      = 200; % Inner points (without boundary), equidistant grid
bounds = [-1, 1];
h      = (bounds(2) - bounds(1)) / n; 

% Init with identity matrix to compensate already for prescribed dirichlet bound. values
A = eye( (n + 1)^2); 
% Equidistant cartesian discretization, used for both x and y
x = linspace(bounds(1), bounds(2), n+1);

bnd_val = 0; % Prescribed Dirichlet boundary condition
b       = bnd_val * ones((n+1)^2, 1);

% Build matrix
for i= 1:n-1
    for j= i * (n + 1) + 2 : i * (n + 1) + n
        A(j,j-(n+1)) = -1/(h^2); % "Far" left sub-diagonal
        A(j, j-1)    = -1/(h^2); % "Close" left sub-diagonal
        A(j,j)       =  4/(h^2);
        A(j,j+1)     = -1/(h^2); % "Close" right sub-diagonal
        A(j,j+(n+1)) = -1/(h^2); % "Far" right sub-diagonal
    end
end 

% Find roughly center point
if mod(n, 2) == 0
  index_center = ( (n + 1)/2 - 0.5 ) * (n + 1) + ( (n + 1)/2 -0.5 );
else
  index_center = ( n/2 - 0.5 ) * (n + 1) + ( n/2 -0.5 );
end
  
% Try the different proposals
b(index_center) = 1/h^2;
%b(index_center) = 3/(4 * h^2);

S = sparse(A); % Accelerate linear system solve
u = S \ b;

% Make approximation ready for 3D plot
U=zeros(n+1);
for i=1:(n+1) 
    U(i, 1:n + 1) = u( (i - 1) * (n + 1) + 1: i * (n + 1) );
end
fprintf('Maximum value of solution vector is: %f\n', max(u) );

% Fundamental solution (solution = 0 only at corners of the [-1, 1]^2 domain)
[X, Y] = meshgrid(x, x);
U_fund = -1 / (2 * pi) * log(sqrt(X.^2 + Y.^2) );

fig = figure;
fig.Position = [100 100 1000 800];

% 3D plots
subplot(2, 2, 1);
surf(x, x, U);
title('Approximation');
shading interp

subplot(2, 2, 2);
surf(x, x, U_fund);
title('Fundamental Solution');
shading interp

% Contour plots
subplot(2, 2, 3);
contour(x, x, U);
title('Approximation');

subplot(2, 2, 4);
contour(x, x, U_fund);
title('Fundamental Solution');

Although it is not quite visible in the Figure, the maximum value of $u$ is for this choice $1.0023$ for a meshwidth of $\Delta x = 0.01$, so of reasonable order.

As a side note: The known fundamental solution has clearly a singularity at $(0, 0)$ and the first derivative at this point does not exist. Thus, $u(x, y)$ is not continuously differentiable, making the Taylor expansion invalid and thus removing all guarantees of the finite difference method. Typically, you would use the Finite Element method for this kind of problem which give you guarantees with much less requirements on $u$.

For this particular problem, you can prove that $$ u (r) = - \frac{1}{2\pi} \ln (r) $$ solves the problem in a weak sense. Define the domain $\Omega = [0, 1]\times [0, 2\pi]$, i.e. the unit disk around the origin. The weak formulation is obtained by multiplying the problem with a test function $v$ which is in this case also chosen to depend on the radius $r$ only, i.e., $w = w(r)$. Then you integrate over the domain to obtain: $$\int_\Omega -\Delta u(r) w(r) \mathrm d \boldsymbol{v} = \int_\Omega \delta w(r) \mathrm d \boldsymbol{v}$$ The right-hand side is by one of the properties of the Dirac Delta just $w(0)$. For the RHS, use integration by parts / reverse product rule in multiple dimensions: $$-\Delta u(r) w(r) = - \nabla \cdot \big(\nabla u(r)\big) w(r) = - \nabla \cdot \Big(\nabla u(r) w(r) \Big) - \Big(- \nabla u \cdot \nabla w \Big)$$ Thus, $$\int_\Omega -\Delta u(r) w(r) \mathrm d \boldsymbol{v} = \int_\Omega - \nabla \cdot \Big(\nabla u(r) w(r) \Big) + \nabla u \cdot \nabla w \mathrm d \boldsymbol{v} $$ Using the Divergence/Gauss Theorem $$\int_\Omega -\Delta u(r) w(r) \mathrm d \boldsymbol{v} = \int_{\partial\Omega} - \Big(\nabla u(r) w(r) \Big) \cdot \boldsymbol{n} \mathrm d s +\int_\Omega \nabla u \cdot \nabla w \mathrm d \boldsymbol{v} $$ For a complete problem formulatio you need, as mentioned already above, boundary conditions. Using homogeneous boundary conditions translates here into $u(r = 1) = 0$. Since $w$ is for various reasons related to functional analysis chosen to have the same property, we have that $w(r = 1) = 0$. Thus, on the boundary $\partial \Omega$ on the unit disk (which is just the circle with radius $1$) you have that $w = 0$, making the integral over the boundary vanish. Therefore, our weak formulation reduced to $$\int_\Omega \nabla u \cdot \nabla w \mathrm d \boldsymbol{v} = w(0) .$$ Due to the radial symmetry of $u$ we shall work in polar coordinates. The $\nabla$-operator is then given by $$\nabla = \begin{pmatrix} \partial_r \\ \frac{1}{r}\partial_\phi \end{pmatrix} $$ and thus $$\nabla u =\begin{pmatrix} - \frac{1}{2\pi} \frac{1}{r} \\ \frac{1}{r} \cdot 0 \end{pmatrix} $$ Thus, we can write \begin{align}\int_\Omega \nabla u \cdot \nabla w \mathrm d \boldsymbol{v} &= \int_0^1 \int_0^{2 \pi} - \frac{1}{2\pi} \frac{1}{r} \partial_r w(r) r \mathrm d \phi \mathrm d r \\ &= - \frac{1}{2\pi} \int_0^1 \int_0^{2 \pi} \partial_r w(r) \mathrm d \phi \mathrm d r \\ &= - \frac{1}{2\pi} \int_0^1 2 \pi \partial_r w(r) \mathrm d r \\ &= - w(r) \bigg\vert_0^1 = - \underbrace{w(r = 1)}_{= 0} - \big(-w(r = 0)\big) = w(0) \checkmark \end{align}