How to solve a second order partial differential equation involving a delta Dirac function? In a mathematical physical problem, I came across the following partial differential equation involving a delta Dirac function:
$$
a \, \frac{\partial^2 w}{\partial x^2}
+ b \, \frac{\partial^2 w}{\partial y^2}
+ \delta^2(x,y) = 0 \, , 
$$
subject to the boundary conditions $w(x = \pm 1, y) = w(x, y = \pm 1) = 0$.
Here $a, b \in \mathbb{R}_+$ and $\delta^2(x,y) = \delta(x)\delta(y)$ is the two-dimensional delta Dirac function.
While solutions for ODEs with delta Dirac functions can readily be obtained using the standard approach, I am not aware of any resolution recipe for PDEs with delta Dirac functions.
Any help or hint is highly desirable and appreciated.
Thank you 
 A: Use an ansatz of the form
$$ w(x,y) = \sum_{n,m=1}^\infty c_{n,m} \sin \left(n\pi \frac{1+x}{2}\right)\sin\left(m\pi \frac{1+y}{2}\right) $$
Decomposing the delta function into its Fourier series gives
$$ \delta(x,y) = \sum_{n,m=1}^\infty \sin \left(\frac{n\pi}{2}\right)\sin\left( \frac{m\pi}{2}\right)\sin \left(n\pi \frac{1+x}{2}\right)\sin\left(m\pi \frac{1+y}{2}\right) $$
Plugging the above expressions into the equation gives
$$ -\left[a\left(\frac{n\pi}{2}\right)^2 + b\left(\frac{m\pi}{2}\right)^2\right]c_{n,m} = -\sin \left(\frac{n\pi}{2}\right)\sin\left( \frac{m\pi}{2}\right) $$
A: *

*The method of images with Dirichlet boundary conditions for a square region $[-1,1]^2$ implies that the 2D Dirac delta distribution $$\delta(x)\delta(y)$$ is part of an alternating Dirac comb/Shah function $$A(x)A(y),$$ where
$$\begin{align}A(x)&~=~\sum_{n\in\mathbb{Z}}(-1)^n\delta(x\!-\!2n)~=~III_4(x)-III_4(x+2)
~=~\frac{1}{4}\sum_{n\in\mathbb{Z}}e^{i\pi n x/2}
-\frac{1}{4}\sum_{n\in\mathbb{Z}}e^{i\pi n(x+2)/2}\cr
&~=~\frac{1}{2}\sum_{n\in\mathbb{Z}}\frac{1-(-1)^n}{2}e^{i\pi n x/2}
~\stackrel{n=2k-1}=~\frac{1}{2}\sum_{k\in\mathbb{Z}}e^{i\pi(k-1/2) x}
~=~\sum_{k\in\mathbb{N}}\cos(\pi(k\!-\!1/2) x)\cr 
&~=~\sum_{k\in\mathbb{N}-\frac{1}{2}}\cos(\pi k x).\end{align}$$

*Therefore the solution to OP's BVP becomes
$$ w(x,y)~=~\frac{1}{\pi^2}\sum_{n,m\in\mathbb{N}-\frac{1}{2}}\frac{\cos(\pi n x)\cos(\pi m y)}{a n^2+b m^2}.$$
We leave it to the reader to analyze convergence properties of the double sum.
A: *

*Idea: If we rescale the coordinates $$(x,y)~=~(\sqrt{a}X,\sqrt{b}Y),$$ the new problem is a 2D electrostatic problem
$$  \left(\frac{\partial^2 }{\partial X^2}+\frac{\partial^2 }{\partial Y^2}\right)w~=~-\frac{1}{\sqrt{ab}}\delta(X)\delta(Y), $$
$$ w(X = \pm 1/\sqrt{a}, Y) ~=~0~=~w(X, Y = \pm1/\sqrt{b}) ,$$
for a rectangular with Dirichlet boundary conditions. The potential $w$ is expected to diverge logarithmically at the location of the point charge. 

*The solution can then be formally obtained via the method of images
$$w~=~-\frac{1}{4\pi\sqrt{ab}} \sum_{n,m\in\mathbb{Z}}(-1)^{n+m}\ln\left\{(X\!-\!\frac{2n}{\sqrt{a}})^2+(Y\!-\!\frac{2m}{\sqrt{b}})^2\right\}$$
$$~=~-\frac{1}{4\pi\sqrt{ab}} \sum_{n,m\in\mathbb{Z}}(-1)^{n+m}\ln\left\{\frac{(x\!-\!2n)^2}{a}+\frac{(y\!-\!2m)^2}{b}\right\}.$$
Hm. The double sum is divergent as the general term doesn't go to zero as $|n|,|m|\to\infty$, cf. below comment by user Winther. We speculate that it may be possible to group alternating terms together to achieve a conditionally convergent series.
