I am attempting to implement a model outlined in this paper:
General magnetostatic shape–shape interactions
Background
This model allows the calculation of magnetostatic interaction energies between objects of arbitrary shape. I am trying to numerically evaluate the integral shown here:
$$E_m=\dfrac{\overline{K}_d}{4\pi^3}\Re\left[\int d^3\mathbf{k} D_1(\mathbf{k})D_1^*(\mathbf{k})\times (\mathbf{\hat{m}_{1}}\cdot \mathbf{\hat{k}}) (\mathbf{\hat{m}_{2}}\cdot \mathbf{\hat{k}})e^{i\mathbf{k}\cdot\mathbf{\rho}}\right]$$
Part of the computation involves the computation of unit vectors in reciprocal space:
$$\mathbf{\hat{k}}= \dfrac{[k_x,k_y,k_z]}{\sqrt{k_x^2+k_y^2+k_z^2}}= \dfrac{[k_x,k_y,k_z]}{k}$$
Where $k_\alpha$ is a frequency in the $x$, $y$, or $z$ axis.
These vectors are grouped in way that you are left with factors like this in the final expression:
$\dfrac{k_x k_z}{k^2}$ or $\dfrac{k_y k_y}{k^2}$
These create singularities at $k=0$ which must be bypassed to complete the numerical calculation.
Question
When I try to evaluate the limit (naively with L'Hopital's rule), I find as $\mathbf{k}$ tends to zero for $\dfrac{1}{k^2}$ or $\dfrac{k_x k_z}{k^2}$, the limit is zero or undefined depending on case. That doesn't seem right, because the $k=0$ term corresponds to the large, DC value of $D(\mathbf{k})$, and that shouldn't simply be zeroed or ignored. A colleague suggested I read into "numerical inverse laplacians" but I wasn't able to find anything that seemed relevant. How can this singularity be bypassed?
Any help is greatly appreciated.