Recall that if $\varphi$ is a scalar field and $\vec{G}$ is a vector field, then: $$\nabla \cdot (\varphi \vec{G}) = \varphi \,\nabla \cdot \vec{G} + (\nabla \varphi) \cdot \vec{G}.$$Make $\varphi = 1/\|\vec{r}\|^3 $ and $\vec{G} = \vec{r}$. Writing $$\varphi(x_1,x_2,x_3) = (x_1^2+x_2^2+x_3^2)^{-3/2},$$we get $$\partial_i\varphi(x_1,x_2,x_3) = -3x_i(x_1^2+x_2^2+x_3^2)^{-5/2} \implies \nabla \varphi = -3\frac{\vec{r}}{\|\vec{r}\|^5}$$So: $$\nabla \cdot\vec{F} = \frac{1}{\|\vec{r}\|^3}3 - 3\frac{1}{\|\vec{r}\|^5}\vec{r}\cdot\vec{r} = \frac{3}{\|\vec{r}\|^3}-\frac{3}{\|\vec{r}\|^3} = 0,$$painlessly.
1) You are right. You can compute the extrema of this function and see how the sources and sinks are distributed and what are their magnitudes.
2) This "stuff" has a name only when it does represents something. I mean, if you are talking about fluid mechanics, then the divergence of the velocity field means compressibility. Since the divergence of a vector field is defined in a point (this is only a mathematical artifact, "because matter is empty") I should have said mass creation or destruction.
To illustrate this, the net amount $\Delta q$ of some physical quantity that is exchanged through any closed surface $S$ within a vector field $\vec{F}$ with constant divergence $\mathrm{div}{\vec{F}}=c$ is proportional only to the volume $V$ this surface encloses:
$$\Delta q =c\,V$$
3) In fact the divergence of a vector field has a definition:
$$\mathrm{div}\vec{F}=\lim_{V\to 0}\frac{1}{V}\int_{\partial V}{\vec{F}\cdot\vec{n}\,d\sigma}$$
where it represents the net exchange of some quantity through a surface enclosing a volume $V$ when it shrinks to $0$.
Of course you can do it that way. You can define a function that tells you the normalised divergence you have if you consider this offers more valuable information to you. Regarding your intervals, yes, you are correct.
4)I am not used to this terminology
5) There are many points of view regarding the answer to this question but mine is that the infinity is not a physical quantity. Infinity means big enough wrt. what is considered (this is related to what I have said before: "matter is empty"). Singularities appear basically when mathematical models cannot describe what is going on in there. For example when assuming a potential solution for the flow past a wing, a singularity appears in the trailing edge, because the viscosity has been obviated (which is a physically incorrect assumption).
Best Answer
First, the divergence in spherical coordinates, expressed in terms of derivatives, would take the form
$$\nabla \cdot \vec A = \frac{1}{r^2} \frac{\partial}{\partial r} [r^2 A^r] + \ldots$$
where $A^r$ is the radial component of the vector field $\vec A$. In this case, that's $1/r^2$, so we naively get 0 for this contribution.
Of course, this formula cannot be valid at the origin (there's a $1/r^2$ in it), so we resort to an alternative definition of divergence, one that is regular at the origin:
$$\nabla \cdot \vec A |_{\vec 0} = \lim_{\delta r \to 0} \frac{1}{4 \pi (\delta r)^3/3} \int_0^{2\pi} \int_0^\pi \vec A \cdot \hat r \, (\delta r)^2 \, \sin \theta \, d\theta \, d\phi$$
This basically follows from the divergence theorem and is a general alternative definition for divergence that eschews derivatives. It also can be done with curl and gradient instead. The limit process, of course, makes this entirely equivalent to a derivative.
Evaluating this definition for $\vec A = \hat r/r^2$ at $r = \delta r$ (as we're integrating over a sphere of this radius, and then taking the limit) gives
$$\left. \nabla \cdot \left( \frac{\hat r}{r^2} \right) \right|_{\vec 0} = \lim_{\delta r \to 0} \frac{3}{\delta r} $$
which obviously diverges (positively, as negative radii have no meaning).
So this is a vector field whose divergence is zero everywhere except the origin, where its divergence...well, diverges. That all certainly sounds like a delta function.
Typically, one uses the divergence theorem directly to verify the stated condition of the delta function: that its integral over any region containing zero is 1. That is, we do
$$4\pi \int_0^r \nabla \cdot \frac{\hat r}{r^2} r^2 \, dr= \int_0^{2\pi} \int_0^\pi \frac{\hat r}{r^2} \cdot \hat r r^2 \sin \theta \, d\theta \, d\phi$$
If $\nabla \cdot \hat r/r^2 = 4 \pi \delta$, then the surface integral on the left should evaluate to $4\pi$ (which it does, obviously). This is not 100% formally sound, but you can always check this by integrating against a test function (a scalar field would do), just as you would in 1d.