Is my interpretation of three dimensional improper integral correct? In Physics/Electrostatics textbook, I am in a situation where we have to find the electric field at a point inside the volume charge distribution. In Cartesian coordinates, we can't do it the usual way because of the integrand singularity. So we use the three dimensional improper integral.
$$\mathbf{E}=\lim\limits_{\epsilon\to 0} \int_{V'-\delta_{\epsilon}}
\rho'\ \dfrac{\mathbf{r}-\mathbf{r'}}{|\mathbf{r}-\mathbf{r'}|^3} dV' \tag1$$
where:
$\mathbf{r'}=(x',y',z')$ is coordinates of source points
$\mathbf{r}=(x,y,z)$ is coordinates of field points
$V'$ is the volume occupied by the charge
$\delta_{\epsilon}$ is an arbitrary volume contained in $V'$ around the singular point $\mathbf{r}=\mathbf{r'}$ with $\epsilon$ being its greatest chord.
$\rho'$ is the charge density and is continuous throughout the volume $V'-\delta_{\epsilon}$
While taking the limit the shape of $\delta_{\epsilon}$ is kept unaltered
From equation $(1)$, we can get the $x$-component of $\mathbf{E}$:
$$E_x=\lim\limits_{\epsilon\to 0} \int_{V'-\delta_{\epsilon}}
\rho'\ \dfrac{x-x'}{|\mathbf{r}-\mathbf{r'}|^3} dV' \tag2$$

I view the steps of solving $E_x$ as follows:


*

*Make a $\delta_{\epsilon}$ cavity (with $\epsilon=a$) contained in $V'$ around the singular point $\mathbf{r}=\mathbf{r'}$. Then take the Riemann integral over $V'-\delta_{\epsilon}$

*Find the function which relates "Riemann integral over $V'-\delta_{\epsilon}$" and "$\epsilon$" over the interval $(0,a]$. For the sake of clarity we can also make a graph of "$\epsilon$" ($x$-axis) and "Riemann integral over $V'-\delta_{\epsilon}$" ($y$-axis) over the interval $(0,a]$

*Find $l$ $\ni$
$\forall \varepsilon > 0, \exists \delta \ni \text{when} |x-0|<\delta, |y-l|< \varepsilon$
Thus $l$ is the solution for $E_x$ in equation $(2)$
I think this is the correct interpretation of equation $(2)$ (as I confirmed with one of my teachers)
I know steps $(1)$ and $(3)$ can be done. Please explain a way to execute step $(2)$.
 A: Your interpretation is correct, so let's talk about practical calculation. If the limit exists, then choosing any cavity shape would get us the correct value, so why not pick a nice one? In practice, one would choose a good $\delta_\epsilon$, say a ball of radius $\epsilon$, to evaluate the integral. 
The ball is a good choice because even if the limit doesn't exist, taking the value in this symmetric way allows us to take the principal value of such an integral, if necessary.
There are two techniques to evaluate the improper integral. For the first, we'll need the following formula: for scalar $f$ and vector $\mathbf{G}$ on $\Omega \subset\subset \mathbb{R}^n$ we have
$$ \int_\Omega f\nabla\cdot \mathbf{G} dr = \int_{\partial \Omega} f\mathbf{G}\cdot \mathbf{\nu} d\sigma - \int_\Omega \nabla f \cdot \mathbf{G} dr $$
where $\mathbf{\nu}$ is the outward normal of the boundary and $d\sigma$ is the surface measure of the boundary. This is the integration by parts formula, and it can be proven with the Divergence theorem.
As an example, take the integral above and let the integrand be $f$ and let $\mathbf{G} = \frac{1}{3}(\mathbf{r}'-\mathbf{r})$. Assuming $\rho\in C^1(\mathbb{R}^3)$, we have
$$\int_{V'-B_\epsilon(\mathbf{r})} \rho(\mathbf{r}')\frac{x-x'}{|\mathbf{r}-\mathbf{r}'|^3} d^3\mathbf{r}' = \int_{\partial V'}\frac{\rho}{3}\frac{x-x'}{|\mathbf{r}-\mathbf{r}'|^3}(\mathbf{r}'-\mathbf{r})\cdot\nu d\sigma$$ $$ - \int_{\partial B_\epsilon(\mathbf{r})}\frac{\rho}{3}\frac{x-x'}{|\mathbf{r}-\mathbf{r}'|^3}(\mathbf{r}'-\mathbf{r})\cdot\frac{\mathbf{r}'-\mathbf{r}}{|\mathbf{r}'-\mathbf{r}|} d\sigma + \int_{V'-B_\epsilon(\mathbf{r})} \nabla\left(\rho(\mathbf{r}')\frac{x-x'}{|\mathbf{r}-\mathbf{r}'|^3}\right)\cdot \frac{\mathbf{r}-\mathbf{r}'}{3}d^3\mathbf{r}' $$
$$ = - \int_{S^2}\frac{\rho}{3}\frac{\epsilon \sin\theta\cos\phi}{\epsilon^2}\epsilon^2 d\Omega + \int_{V'-B_\epsilon(\mathbf{r})} \nabla\left(\rho(\mathbf{r}')\frac{x-x'}{|\mathbf{r}-\mathbf{r}'|^3}\right)\cdot \frac{\mathbf{r}-\mathbf{r}'}{3}d^3\mathbf{r}' $$
where we used the fact that the charge distribution vanishes at the boundary, converted to spherical coordinates centered at $\mathbf{r}$ for the remaining surface integral, and $d\Omega$ is angular measure on the unit sphere. Notice that as $\epsilon\to 0$, that boundary integral vanishes (due to $\rho$ being continuous) and we are only left with the volume integral on the right. 
Ordinarily this is not the way I recommend this strategy be implemented $-$ try to find a vector field whose divergence contains the singularity, rather than trying to continue to differentiate it (unless the derivative eventually vanishes, like for the Laplacian of $|\mathbf{r}-\mathbf{r}'|^{-1}$). But computational tricks depend on the exact form of $\rho$. The point of this was to give an example of how these integration by parts computations tend to go, and where to take the limit.
The other technique is to move to a coordinate system where the singularity vanishes. Take the original integral and cut it up like this: 
$$\int_{V'-B_\epsilon(\mathbf{r})} \rho(\mathbf{r}')\frac{x-x'}{|\mathbf{r}-\mathbf{r}'|^3} d^3\mathbf{r}' = \int_{V'-B_R(\mathbf{r})} \rho(\mathbf{r}')\frac{x-x'}{|\mathbf{r}-\mathbf{r}'|^3} d^3\mathbf{r}' + \int_{B_R(\mathbf{r})-B_\epsilon(\mathbf{r})} \rho(\mathbf{r}')\frac{x-x'}{|\mathbf{r}-\mathbf{r}'|^3} d^3\mathbf{r}'$$
The first integral is not improper so we'll ignore it. Now convert the second integral into spherical coordinates centered at $\mathbf{r}$:
$$\int_0^{2\pi}\int_0^\pi\int_\epsilon^R \rho(\mathbf{r}')\frac{-r''\sin\theta\cos\phi}{r''^3} r''^2\sin\theta dr''d\theta d\phi \to \int_0^{2\pi}\int_0^\pi\int_0^R -\rho(\mathbf{r}')\sin^2\theta\cos\phi dr''d\theta d\phi$$
where, since $\rho$ is continuous, there is now no singularity at $r''=0$ because of the Jacobian. 
Hopefully that clears up the computational side of improper integrals in higher dimensions.
