Divergence in spherical coordinates vs. cartesian coordinates The divergence in spherical coordinates is derived as
$$\nabla\vec A=\frac{1}{r^2}\frac{\partial(r^2A_r)}{\partial r}+\frac{1}{r\sin\theta}\frac{\partial}{\partial\theta}(A_\theta\sin\theta)+\frac{1}{r\sin\theta}\frac{\partial A_\varphi}{\partial\varphi}$$
However, this seems to fail at a simple example. Let's take the identity $\vec F(\vec r):=\vec r$ for example. The divergence of the vector field $\vec F(\vec r)=\begin{pmatrix}x\\y\\z\end{pmatrix}$ in cartesian coordinates is
$$\nabla\vec F(\vec r)=\frac{\partial}{\partial x}F_x(\vec r)+\frac{\partial}{\partial y}F_y(\vec r)+\frac{\partial}{\partial z}F_z(\vec r)=\frac{\partial}{\partial x}x+\frac{\partial}{\partial y}y+\frac{\partial}{\partial z}z=3$$
But in spherical coordinates the vector field would be fined as $\vec F(\vec r)=\begin{pmatrix}r\\\theta\\\varphi\end{pmatrix}$, so the divergence would be:
$$\nabla\vec F(\vec r)=\frac{1}{r^2}\frac{\partial}{\partial r}(r^2\cdot r)+\frac{1}{r\sin\theta}\frac{\partial}{\partial\theta}(\theta\sin\theta)+\frac{1}{r\sin\theta}\frac{\partial}{\partial\varphi}(\varphi)=3+\frac{\sin\theta+\theta\cos\theta}{r\sin\theta}+\frac{1}{r\sin\theta}$$
Which certainly isn't 3:

(Graphic from WolframAlpha)
So what am I doing wrong here?

Edit I successfully managed to completely confuse myself. I think I misunderstood the meaning of $A_r, A_\theta, A_\varphi$ - could anybody explain it to me? In my understanding the following should be true:
$$A_r=\vec A\cdot\hat e_r\qquad A_\theta=\vec A\cdot\hat e_\theta\qquad A_\varphi=\vec A\cdot\hat e_\varphi$$
With the unit vectors for spherical coordinates being as follows:
$$\hat {e} _{r}={\begin{pmatrix}\sin \theta \cos \varphi \\\sin \theta \sin \varphi \\\cos \theta \end{pmatrix}},\qquad \hat {e} _{\theta }={\begin{pmatrix}\cos \theta \cos \varphi \\\cos \theta \sin \varphi \\-\sin \theta \end{pmatrix}},\qquad \hat {e} _{\varphi }={\begin{pmatrix}-\sin \varphi \\\cos \varphi \\0\end{pmatrix}}$$
However, this doesn't make much more sense, as - if I understood it correctly, $\varphi$ and $\theta$ would be $45^\circ$ (and $r=1$) which results in weird basis vectors and doesn't seem to help me.
 A: I am updating this answer to try to address the edited version of the question.
A nice thing about the conventional $(x,y,z)$ Cartesian coordinates is everything works the same way. In fact, everything works so much the same way using the same three coordinates in the same way all the time in Cartesian coordinates--points in space, vectors between points, field vectors--that it may be surprising that in just about any other coordinate system different things sometimes work differently from each other.
In Cartesian coordinates, you can identify the components of a point's coordinates with the components of the point's position vector.
And you can get the vector sum of two of those vectors by adding the coordinates:
$$
 \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}
 + \begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix}
 = \begin{pmatrix} 1+0 \\ 0+1 \\ 0+0 \end{pmatrix}
 = \begin{pmatrix} 1 \\ 1 \\ 0 \end{pmatrix}.
$$
But if you try to describe a vectors by treating them as position vectors and using the spherical coordinates of the points whose positions are given by the vectors, the left side of the equation above becomes
$$
 \begin{pmatrix} 1 \\ \pi/2 \\ 0 \end{pmatrix}
 + \begin{pmatrix} 1 \\ \pi/2 \\ \pi/2 \end{pmatrix},
$$
while the right-hand side of the equation becomes
$$
 \begin{pmatrix} \sqrt2 \\ \pi/2 \\ \pi/4 \end{pmatrix}.
$$
So spherical coordinates are really awkward to do vector operations in. You just can't do it like you can with Cartesian coordinates.
Another thing that's going on with these formulas is we have to make sense of things like $\frac\partial{\partial\theta}.$
If you have something that can vary as $r,$ $\theta,$ and $\varphi$ vary, then
$\frac\partial{\partial\theta}$ is looking for what happens when we hold $r$ and $\phi$ constant but vary $\theta.$
What happens to the location defined by the coordinates $(r, \theta, \varphi)$
if we hold $r$ and $\varphi$ constant but increase $\theta$ is that we move in a "downward" direction perpendicular to the position vector.
Given a starting point $(r, \theta, \varphi)$ for this motion, we define
$\hat{\boldsymbol \theta}$ as the unit vector in that direction.
Similarly, there is a unit vector $\hat{\boldsymbol \varphi}$ in the direction we would travel if we hold $r$ and $\theta$ constant and increase $\varphi.$
There is also a unit vector $\hat{\mathbf r}$ in the direction of constant $\theta$ and $\varphi$ but increasing $r.$
I name these vectors $\hat{\mathbf r},$ $\hat{\boldsymbol \theta},$ and
$\hat{\boldsymbol \varphi}$ because that's how they're named on the first Wikipedia page the question linked to. These are the same vectors named
$\hat e_r,$ $\hat e_\theta,$ and $\hat e_\varphi$ elsewhere.
Unlike the unit vectors $\hat{\mathbf x},$ $\hat{\mathbf y},$ and
$\hat{\mathbf z},$ which are always parallel to the $x,$ $y,$ and $z$ axes,
the unit vectors $\hat{\mathbf r},$ $\hat{\boldsymbol \theta},$ and
$\hat{\boldsymbol \varphi}$ depend on the coordinates of a point 
$(r, \theta, \varphi)$ relative to which those vectors are defined.
Specifically, 
$$\begin{align}
\hat{\mathbf r} &= \sin\theta \cos\varphi \hat{\mathbf x}
 + \sin\theta\sin\varphi \hat{\mathbf y} + \cos\theta \hat{\mathbf z}, \\
\hat{\boldsymbol \theta} &= \cos\theta \cos\varphi \hat{\mathbf x}
 + \cos\theta \sin\varphi \hat{\mathbf y} - \sin\theta \hat{\mathbf z}, \\
\hat{\boldsymbol \varphi} &= -\sin\varphi \hat{\mathbf x} + \cos\varphi \hat{\mathbf y}.
\end{align}
$$
You can visualize a physical model for these vectors by supposing that the Earth is spherical and is centered at the origin of the coordinate system. Then at any point on the surface of the Earth, $r$ is the Earth's radius,
$\theta$ is the co-latitude (like latitude, but it goes from zero at the north pole to $180$ degrees at the south pole), and $\varphi$ is the longitude.
Then $\hat{\boldsymbol \theta}$ is in the direction of increasing co-latitude, so it points due south;
$\hat{\boldsymbol \varphi}$ is in the direction of increasing longitude, so it points due east;
and $\hat{\mathbf r}$ is in the direction of increasing radius, so it points directly away from the center.
In particular, $\hat{\mathbf r}$ is a unit vector in the exact same direction as the position vector $\vec r,$
and the other two unit vectors are orthogonal to $\hat{\mathbf r}$ and to each other.
This clearly does not work as the basis of a set of position vectors for all the points in a three-dimensional space. If you wanted the vector describing the direction and distance from point $P$ to point $Q,$ do you use the vectors 
$\hat{\mathbf r},$ $\hat{\boldsymbol \theta},$ and
$\hat{\boldsymbol \varphi}$ according to how they are defined relative to $P$ or according to how they are defined relative to $Q$?
The answer is (usually) neither, you don't typically use those vectors for that kind of calculation at all.
What those vectors are good for is as an orthonormal basis for things that occur locally at a point: a vector field at that point, the gradient of that field, the divergence of that field, and so forth.
Now if you have a vector field with the value $\vec A$ at some point with spherical coordinates 
$(r, \theta, \varphi),$ then we can break that vector down into orthogonal components exactly as you do:
\begin{align}
A_r &= \vec A\cdot\hat{\mathbf r}, \\
A_\theta &= \vec A\cdot\hat{\boldsymbol \theta}, \\
A_\varphi &= \vec A\cdot\hat{\boldsymbol \varphi}.
\end{align}
Now consider the case where $\vec A = \vec r.$ Then $\vec A$ is in the exact same direction as
$\hat{\mathbf r},$ and therefore
$$ A_r = \vec A\cdot\hat{\mathbf r} = \lVert \vec A\rVert \lVert \hat{\mathbf r}\rVert = r. $$
On the other hand, $\vec A$ is orthogonal to each of the vectors
$\hat{\boldsymbol \theta}$ and $\hat{\boldsymbol \varphi},$ so
$$ \vec A\cdot\hat{\boldsymbol \theta} = \vec A\cdot\hat{\boldsymbol \varphi} = 0.$$
In summary, $A_r = r$ and $A_\theta = A_\varphi = 0.$
Moreover, no matter which way you go, $A_\theta$ and $A_\varphi$ both remain zero, so 
$$
\frac{\partial}{\partial\theta}(A_\theta\sin\theta)
 = \frac{\partial A_\varphi}{\partial\varphi} = 0.
$$
If you are measuring the field at the point $(1,1,1),$ then $\hat{\mathbf r}$ at that point is 
$\hat{\mathbf r} = \frac1{\sqrt3}\hat{\mathbf x} + 
\frac1{\sqrt3}\hat{\mathbf y} + \frac1{\sqrt3}\hat{\mathbf y},$
and therefore 
$\sqrt3\hat{\mathbf r} = \hat{\mathbf x} + \hat{\mathbf y} + \hat{\mathbf z},$ or in the usual $(x,y,z)$ Cartesian system, the vector $(1,1,1).$
Plugging in the correct values for all three components, the divergence formula gives the result $3.$
Now if you had a field that was not simply a radial field always going directly outward from the center or directly inward toward the center, that is, if there were some field vectors making other angles with the radial vectors, then you would see some non-zero values of $A_\theta$ and $A_\phi,$ and there would be a reason you needed the rest of the formula for divergence. You just happen to have chosen an example in which you have eliminated those components.
If you want a more interesting example, you might try a vector field such as
$\vec A = x \hat{\mathbf x},$ which has non-zero derivatives in all three directions at most points.
A: It is sadly common when dealing with different coordinate systems, such as in this case, not to distinguish the function under consideration from its pullback under a coordinate transformation.
Pullback.
Let me introduce some terminology. Suppose we have sets $A$, $B$, and $C$ along with maps
$$A \xrightarrow{\Phi} B \xrightarrow{u} C.$$
We may then form the composition $u \circ \Phi$ of $u$ with $\Phi$, which is the mapping of $A$ into $C$,
$$A \xrightarrow{u \circ \Phi} C,$$
defined by letting
$$(u \circ \Phi)(a) = u(\Phi(a)) \quad\text{for all}\quad a \in A.$$
We call $u \circ \Phi$ the pullback of $u$ under $\Phi$, sometimes denoted $\Phi^\ast u$. Below we will also denote this map simply by $\tilde u$, thus
$$\tilde u = \Phi^\ast u = u \circ \Phi$$
will all denote the same map in the below considerations. (The trouble arises when people do not distinguish between $u$ and $\tilde u$.)
Spherical Coordinate Example.
Consider now the subset $A$ of $\mathbf{R}^3$ given by
$$A = [0, \infty) \times [0, \pi] \times [0, 2\pi] = \bigl\{ (r, \theta, \phi) : r \ge 0, 0 \le \theta \le \pi, 0 \le \phi \le 2\pi \bigr\}$$
along with the mapping $\Phi$ of $A$ into $\mathbf{R}^3$ given by
$$\Phi(r, \theta, \phi) = \begin{bmatrix} r\sin\theta\cos\phi \\ r\sin\theta\sin\phi \\ r\cos\theta \end{bmatrix}.$$
If for instance $u$ is a mapping of $\mathbf{R}^3$ into $\mathbf{R}$, we call the pullback $\tilde u = u \circ \Phi$ of $u$ under $\Phi$ the representation of $u$ in spherical coordinates.
As an example, suppose $u$ is defined by letting
$$u(x, y, z) = (x^2 + y^2 + z^2)^{\frac 12} \quad\text{for all}\quad (x, y, z) \in \mathbf{R}^3,$$
and let us compute the representation of $u$ in spherical coordinates. For a point $(r, \theta, \phi) \in A$, one has
$$\begin{align*}
\tilde u(r, \theta, \phi)
&= u(\Phi(r, \theta, \phi))
= u\bigl( r\sin\theta\cos\phi, r\sin\theta\sin\phi, r\cos\theta \bigr) \\
&= \bigl( (r\sin\theta\cos\phi)^2 + (r\sin\theta\sin\phi)^2 + (r\cos\theta)^2 \bigr)^{\frac 12} \\
&= r.
\end{align*}$$
Thus, $u$ is the mapping on $\mathbf{R}^3$ sending the point $(x,y,z)$ to $(x^2 + y^2 + z^2)^{1/2}$, whereas $\tilde u$ is the mapping on $A$ sending the point $(r, \theta, \phi)$ to $r$, i.e., projection on the first factor.
Observe that $u$ looks a lot simpler in spherical coordinates! It is therefore of interest to see if we can compute such quantities as the gradient $\nabla u$ of $u$ in terms of the gradient $\nabla\tilde u$ of $\tilde u$ and the mapping $\Phi$.
Vector Fields in Spherical Coordinates.
The representation of a vector field in spherical coordinates is worked out in the same way. Consider your example $F$ mapping a point $(x,y,z)$ in $\mathbf{R}^3$ to itself:
$$F(x,y,z) = \begin{bmatrix} x \\ y \\ z \end{bmatrix}.$$
Its pullback under $\Phi$ at a point $(r, \theta, \phi) \in A$ is given by
$$\begin{align*}
\tilde F(r, \theta, \phi)
&= F(\Phi(r, \theta, \phi)) \\
&= F\bigl( r\sin\theta\cos\phi, r\sin\theta\sin\phi, r\cos\theta \bigr) \\
&= \begin{bmatrix} r\sin\theta\cos\phi \\ r\sin\theta\sin\phi \\ r\cos\theta \end{bmatrix}.
\end{align*}$$
Note carefully that although $F$ is the identity mapping on $\mathbf{R}^3$, the representation $\tilde F$ of $F$ in spherical coordinates is not the identity mapping on $A$! However, an interesting question is whether one can recover $\operatorname{div} F$ from $\tilde F$ and $\Phi$, which is what the formula you provided is meant to do.
As for the notation «$F_r$, $F_\theta$, $F_\phi$», one means the projection of $\tilde F$ along the directions of the unit vectors $e_r$, $e_\theta$, and $e_\phi$, respectively. (It would be better to denote these by $\tilde F_r$, $\tilde F_\theta$, $\tilde F_\phi$, but we will not do so below.) Thus at a point $(r, \theta, \phi)$, one means for instance
$$F_r(r, \theta, \phi) = \tilde F(r, \theta, \phi) \cdot e_r(r, \theta, \phi),$$
where we have explicitly evaluated $e_r$ at $(r, \theta, \phi)$.
It is an exercise to show that
$$F_r(r, \theta, \phi) = r, \quad F_\theta(r, \theta, \phi) = 0, \quad F_\phi(r, \theta, \phi) = 0.$$
Using this, the formula you provide for the divergence of $F$ in spherical coordinates becomes
$$\begin{align*}
(\operatorname{div} F)(r\sin\theta\cos\phi, r\sin\theta\sin\phi, r\cos\theta)
&= (\operatorname{div} F)(\Phi(r, \theta, \phi)) \\
&= \frac{1}{r^2}\frac{\partial}{\partial r} \bigl( r^2 F_r(r, \theta, \phi) \bigr) \\
&= \frac{1}{r^2} \frac{\partial}{\partial r}\bigl( r^3 \bigr) \\
&= 3,
\end{align*}$$
which is correct.
A: You didn't transform your coordinates properly:
$$
\left [\begin{array}{c}
x\\
y\\
z\\
\end{array} \right]  \longrightarrow  \left[\begin{array}{c}
\rho \sin\theta\cos\phi\\
\rho\sin\theta\sin\phi\\
\rho\cos\theta\\
\end{array}\right]
$$
