I was reading an answer to a stack exchange post titled Is the electric field of a volume charge distribution well defined? . That answer is shown in the image below:
Now I make a comparison between $\displaystyle\iiint_V\frac{\rho (\vec{r'})(\vec{r}-\vec{r'})}{\lvert (\vec{r}-\vec{r'})^{3}\rvert}d^{3}\vec{r} \text{ and} \int^2_{-1} \dfrac{1}{x^4}dx$. They both are similar in the sense that both contain a singularity point.
First let us consider the former:
Here we consider an infinitesimal volume $dV$ containing the singularity point.
$\displaystyle\iiint_{dV} \frac{\rho (\vec{r'})(\vec{r}-\vec{r'})}{\lvert (\vec{r}-\vec{r'})^{3}\rvert}d^{3}\vec{r}$
Following the reasoning presented in that answer, we see that this integral vanishes in the limit as $\vec{r}\rightarrow\vec{r'}$
Let us now consider the latter:
Here we consider an infinitesimal element of $x$-axis containing the singularity point.
$\displaystyle \int^{+dx}_{-dx}\dfrac{1}{x^4}dx=\left[ \dfrac{x^{-3}}{-3} \right]^{+dx}_{-dx}= -\dfrac{1}{3} \left[ \dfrac{1}{x^3} \right]^{+dx}_{-dx}= -\dfrac{1}{3} \left[ \dfrac{1}{dx^3}-\dfrac{1}{(-dx)^3} \right] =- \dfrac{2}{3 \ dx^3}= -\infty$
But a problem arises here:
$\displaystyle \int^2_{-1} \dfrac{1}{x^4}dx= -\dfrac{3}{8}$
How come a portion of our area $\displaystyle \int^{+dx}_{-dx}\dfrac{1}{x^4}dx$ be greater than our full area $\displaystyle \int^2_{-1} \dfrac{1}{x^4}dx$?
What both integrals have in common is that strictly speaking, both are not well defined as they contain a singularity. However there are methods to deal with the singularity (see below). But what is not valid is to just ignore the singularity and use “the” antiderivative at the end points.
The way you can treat such a singularity is to exclude a symmetric interval around it from the integral, and then take the limit of that interval's length to zero. That is, you actually calculate e.g. for your last integral $$ \begin{aligned} \lim_{\epsilon\to0}\left(\int_{-1}^{-\epsilon}\frac{1}{x^4}\,\mathrm dx + \int_{\epsilon}^2\frac{1}{x^4}\,\mathrm dx\right) &=\lim_{\epsilon\to0}\left(\left.\left(-\frac{1}{3(x)^3}\right)\right|_{-1}^\epsilon + \left.\left(-\frac{1}{3(x)^3}\right)\right|_{\epsilon}^1\right)\\ &=\lim_{\epsilon\to0}\left(-\frac{1}{3(-\epsilon)^3}+\frac{1}{3(-1)^3}-\frac{1}{3\cdot2^3}+\frac{1}{3\epsilon^3}\right)\\ &=\lim_{\epsilon\to0}\left(-\frac{3}{8}+\frac{2}{3\epsilon^3}\right) = \infty \end{aligned} $$ And if you look at the graph of the function $\frac{1}{x^4}$ you'll see that it is always positive, so it would be absurd that the integral over it (i.e. the area below it) would be negative.
Now what is different in the integral from the physics.SE answer? Well, let's get rid of $\rho(\mathbb r)$ (since the assumption is that you can treat it as constant), assume that $\mathbb r'=0$, and reduce it to just one dimension, replacing $\mathbb r$ by $x$. Also, since we are integrating over a sphere, we get a symmetric integral. Then we get $$\int_{-a}^a \frac{x}{|x|^3}\,\mathrm dx$$ Note that the function under the integrand is odd, that is, when replacing $x$ by $-x$ it changes sign. Thus if you do the same replacement as before, the left integral is the exact negative of the right integral, so we get exactly zero. Of course an exact zero remains an exact zero when taking the limit of $\epsilon\to0$.
Note that this would be different if the absolute value bars in the denominator had not been there, or if the numerator would have the absolute value bars as well; in that case, the integral would diverge to $+\infty$ as well.