How to solve a second order partial differential equation involving a delta Dirac function?

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

1. 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}$$

2. 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.