Why is Divergence of a vector field which is decreasing in magnitude as we move away from origin positive at points other than origin? 
divergence of $\frac{\hat{r}}{r}$ is positive, $\frac{\hat{r}}{r^2}$ is zero and $\frac{\hat{r}}{r^3}$ is negative at points other than origin.

As I have studied, divergence in 3D space tells whether there are sources or sinks present.My initial thinking was that since the magnitude of the vectors are decreasing, so the divergence should be negative for all the three cases but that's not the case.It is not like more vector field is being generated, the same number of vector field passes through the spherical cross section of radius $r+dr$ that is being passed through the cross section of radius $r$.
My Question is can anyone provide the physical intuition of why the divergence is positive for the first vector field, zero for the second and negative for the third?
 A: First let build up some intuition by looking at the divergence in 1D. Imagine you have a river which has a very uniform flow so you can approximate the flow as 1 dimensionsal. The flow is also static in time. This particular river flows in the positive x direction so $v(x)$ is positive and assume at a particular location we have that $\frac{\partial v}{\partial x}>0$. This means the divergence can be approximated by
$$\nabla\cdot\vec v=\frac{\partial v}{\partial x}\approx\frac{1}{\Delta x}(v(x+\Delta x)-v(x))>0$$
We can conclude that $v(x+\Delta x)> v(x)$. The quantity $v(x)\times\text{cross section}$ gives how much volume is passing through the plane at $x$ at each second. We must conclude that more volume is passing through $v(x+\Delta x)$ than through $v(x)$. (we assume constant cross section). The density of the water in the region $[x,x+\Delta x]$ must be decreasing as the water flows through this section of river given that no water is destroyed/created. So if you imagine a vector field $\vec v$ as being the flow of water/particles then $\nabla\cdot \vec v>0$ means the density decreases locally and conversely when $\nabla\cdot \vec v<0$ the density increases locally.
So now that we have built up some intuition I want to show it visually. It's better to draw this in 2D so let's consider the divergence of these functions in 2D:
\begin{align}
\nabla\cdot\left(\hat r r\right)&=2>0\\
\nabla\cdot\left(\frac{\hat r}{r}\right)&=0\\
\nabla\cdot\left(\frac{\hat r}{r^3}\right)&=-\frac{2}{r^4}<0
\end{align}
I picked these functions to give the best visual results.
In the following animations you see particles following a vector field as velocity. The vector fields are, from top to bottom, given by $\vec v(\vec r)\propto\vec r, \frac{\hat r}{r}, \frac{\hat r}{r^3}$. They are also scaled a bit because otherwise some of these animations would be too quick or too slow.



It may be a bit hard to see but you can tell that in the top picture (positive divergence) the density goes down as the particles move outwards. This effect is a combination of the particles becoming more spread out and the particles speeding up, both of which increase divergence. In the bottom picture you see the particles bunching up. The effect of slowing down is so strong that it completely overpowers the spreading out effect. In the middle picture the spreading out and slowing down are perfectly balanced which results in a nice uniform density.
These effects can be translated to 3D but the particular dimensions change, the effect of spreading out is stronger in higher dimensions.
A: The problem with the divergence of the fields you wrote is that it is ill-defined in the origin. So whatever you find is valid only if $r\neq0$ and we need to manually add the value of the divergence in the origin. How do we do it? We use the fact that the integral of the divergence in the volume is equal to the flux of the vector field on a surface which encloses that volume.
We start by computing the flux $\Phi$ of your vector fields on a spherical shell $S$ of radius $R$ i.e
$$\Phi = \int dS {1\over R^n}$$
where I used the fact that the vector field is always perpendicular to the surface so we can just integrate its value at $r=R$ on the surface $S$. Of course, in spherical coordinates, $dS=R^2\sin(\theta)d\theta d\phi$ hence
$$\Phi = R^{2-n}\int  \sin(\theta)d\theta d\phi = 4\pi R^{2-n} $$
where the integral I did is just the solid angle $4\pi$.
This must be correspond to the integral of the divergence inside the volume.
As you can see, the flux on the surface is not always the same and can depend on $R$. This is because, except the $n=2$ case, the other fields decrease too fast / not fast enough for the flux to always be the same at whatever distance (the surface scales as $\sim R^2$ so you need to compensate for that or you will get an $R$-dependency)
For n=2 the flux is $4\pi$. But if we integrate the divergence, we get 0. To make things work out we need to define a Dirac's delta function in the origin that fixes thing. A bit tricky but works. See comments to your questions.
For n=1 the flux is $4\pi R$ and the divergence is $1/r^2$. If we integrate the divergence over the volume we get
$$4\pi\int^R r^2dr\; {1\over r^2} = 4\pi R$$ so all good. The flux (and hence also the integral of the divergence over the volume) are $R$ dependent but both quantities are the same. All good. We don't even have to care about the behavior in the origin, it just works..
For $n=3$ (and $n>2$ in general) things get tricky. We know, as you mentioned, that the flux (and hence the integral of the div) must be positive and indeed we get $\Phi = 4\pi /R>0$ (nut decreasing with $R$ as the field goes to 0 very quick). However if we try to do the integral in the volume we find
$$-4\pi\int ^R r^2dr {1\over r^4}=\lim_{\epsilon\to 0} 4\pi\left({1\over R} -{1\over \epsilon}\right)$$
where I used $\epsilon$ to exclude the origin. This is negative and does not converge. Again, the problem is that we don't know what is happening in the origin. We just know that there must be something so big (so infinite!) that it compensates the infinity of the above integral and also its being negative, it must be something like $\sim +4\pi /\epsilon$ with $\epsilon\to 0$ (I don't actually know, but I assume in the origin there is something like $4\pi\delta(r)/r$)
More in general, the problem is that the vector field is decreasing so fast that the only relevant contribution is in the origin and we don't know it. That will fix the "minus sign" paradox. This results in the divergence pointing towards the origin because in the origin there is something huge that compensates.
Your confusion comes from the fact that you consider the divergence a measure of the direction of the field (they all point radially outwards) where instead is a measure of what behavior the flux has as you move away from the origin (bigger values of $R$): constant if $n=0$, increasing if $n<2$, decreasing if $n>2$ and that is why the divergence is 0 [constant flux], positive [increasing flux] or negative [decreasing flux]. Roughly.
The divergence is also not a measure of how much vector field is generating but rather of how much flux the source in the origin provides.
