Simplify ignores DiracDelta contribution from Coulomb sources

Question

Executing the following line,

Simplify[Div[Grad[-1/Sqrt[{x, y, z}.{x, y, z}], {x, y, z}], {x, y, z}]]

yields 0. However, this is incorrect, since $\nabla\cdot\nabla\left(|\mathbf{r}|^{-1}\right)=-4\pi\delta(\mathbf{r})$ in 3-space.

I understand that this problem is coming from the fact that Simplify assumes $\mathbf{r}\neq\mathbf{0}$. However, I don't know how to coax Mathematica into giving the correct result without essentially having to resort to doing mathematics by hand, which kind of defeats the purpose of using a computer algebra system. In this case the result is obviously wrong so no harm is done, but for more complicated vector expressions I have no idea whether Simplify is going to give wrong answers without warning when applied to functions with spatial singularities.

Is there any automated way to handle such singularities to give correct answers in Mathematica?

Answer 1

Here's one way, which certainly doesn't handle all singularities. First we gather the denominator of every subexpression.

allDens[expr_] := Union[Cases[expr, e_ :> Denominator[e], {0, Infinity}]]

From here, we'll just throw a warning that each denominator is assumed to be nonzero.

acknowledgeSingularitiesSimplify[expr_] := Module[{dens},
  dens = Thread[allDens[expr] == 0];
  dens = Cases[Union[Simplify[dens]], _Equal];

  ConditionalExpression[Simplify[expr], Reduce[Not[Or @@ dens]]]
]

acknowledgeSingularitiesSimplify[
  Div[Grad[-1/Sqrt[{x, y, z}.{x, y, z}], {x, y, z}], {x, y, z}]]

(* ConditionalExpression[0, Sqrt[x^2 + y^2 + z^2] != 0] *)

Another thing you might look into, to fully be integrated with Simplify, is the option TransformationFunction.