Derivation of divergence in spherical coordinates from the divergence theorem I'm trying to find the expression of the divergence of a vector field $\vec{E}$ in spherical coordinates from the theorem :
$$\iint_{S(V)}(\vec{E}.\vec{n})dS = \iiint_{V}div(\vec{E})dV$$
but if I write $\vec{E}$ in spherical coordinates:
$$\vec{E} = E_r\vec{e_r}+E_{\phi}\vec{e_{\phi}}+E_{\theta}\vec{e_{\theta}}$$
and if I consider a spherical volume and its surface, I find that $\vec{n} = \vec{e_r}$ since $\vec{n}$ is orthogonal to the spherical surface at any point... So I'm left with $(\vec{E}.\vec{n})  = E_r$ and 
$$\iint_{S(V)}E_rdS = \iiint_{V}div(\vec{E})dV$$
I don't understand how I'm supposed to get to
 $$div(\vec{E}) = \frac{1}{r^2}\frac{\partial (r^2 E_r)}{\partial r} + \frac{1}{r sin\theta}\frac{\partial E_\phi}{\partial\phi}+\frac{1}{r sin\theta}\frac{\partial(sin\theta E_{\theta})}{\partial \theta}$$
 A: Consider a small volume around a point $(r, \theta, \phi)$. That is, try to compute this:
$$\int_{r}^{r+\delta r} \int_{\theta}^{\theta + \delta \theta} \int_{\phi}^{\phi + \delta \phi} \nabla' \cdot \vec E(\vec r') \, {r'}^2 \sin \theta' \, dr' \, d\theta' \, d\phi'$$
You should be able to argue that the lowest-order term in $\delta r \, \delta \theta \, \delta \phi$ is $[\nabla \cdot \vec E(r, \theta, \phi)] r^2 \sin \theta \, \delta r \, \delta \theta \, \delta \phi$.
Now, construct the corresponding surface integral from the divergence theorem.  Your surface will have six smooth pieces (like a cube, but the surfaces are curved, as they follow the coordinate lines).  Each surface's normal direction is another coordinate direction, and as such, you only need to consider one component on each face (for instance, on the $\theta \phi$ surfaces, you only consider the radial component, as the others must contribute nothing to the dot product).
I'll compute two of the faces directly:
$$\int_{\theta}^{\theta + \delta \theta} \int_{\phi}^{\phi + \delta \phi}  [\vec E(r+\delta r, \theta', \phi')(r+\delta r)^2 - \vec E(r, \theta', \phi')r^2] \cdot \hat r \sin \theta' \, d\theta' \, d\phi'$$
This is a computation for two of the six faces of this not-exactly-cube-shaped surface.  The $r+\delta r$ part corresponds to the face furthest from the origin, and the $r$ part corresponds to the face closest to the origin.
Again, consider the lowest order terms in $\delta r \, \delta \theta \, \delta \phi$.  Write $\vec E(r + \delta r, \theta, \phi) = \vec E(r, \theta, \phi) + \partial \vec E/\partial r |_{r, \theta, \phi} \delta r + O(\delta r^2)$, and we get
$$\left[ 2r E_r(r, \theta, \phi) + r^2 \frac{\partial E_r}{\partial r}|_{r, \theta, \phi} \right] \sin \theta \, \delta r \, \delta \theta \, \delta \phi$$
Observe that $\frac{\partial{(r^2 \vec E)}}{\partial r} = 2r\vec E + r^2 \partial \vec E/\partial r$, and I think you can see where this is going.
Repeat this procedure on the other four faces of the pseudo-cube surface, keeping in mind what the normal vectors are (they're not all $\hat r$!).  If you're uncomfortable thinking about keeping around only these terms that are of order $\delta r \, \delta \theta \, \delta \phi$, you can think of it as a limit instead.  You can also use a symmetric integration volume instead, if you're comfortable keeping around factors of 2.
Indeed, this process is sometimes used as a definition of the divergence and other differential operations.  Doing so makes the fundamental theorem of calculus (divergence theorem, etc.) trivial consequences of the definitions.
A: I assume you have some familiarity with the calculus of differential forms in $\mathbb{R}^3$  $\newcommand{\bR}{\mathbb{R}}$ $\newcommand{\ii}{\boldsymbol{i}}$ $\newcommand{\jj}{\boldsymbol{j}}$ $\newcommand{\kk}{\boldsymbol{k}}$ Suppose 
$$E= E_x\ii+E_y\jj+E_z\kk. $$
We associate to it the differential $2$-form
$$\Phi_E= E_x dy\wedge dz+E_y dz\wedge dx +E_z dx\wedge dy.$$
Convert it to spherical coordinates, then take the exterior derivative. You will  obtain something that has the form $\newcommand{\vfi}{\varphi}$
$$w(r, \theta, \vfi) dr\wedge d\theta\wedge d\vfi. $$
The divergence of $E$  is  the quantity
$$\frac{w(r,\theta,\vfi)}{r^2\sin\vfi}. $$
A: Here's a way of calculating the divergence.
First, some preliminaries. The first thing I'll do is calculate the partial derivative operators $\partial_x,\partial_y,\partial_z$ in terms of $\partial_r, \partial_\theta, \partial_\varphi$. To do this I'll use the chain rule. Take a function $v:\Bbb R^3\to\Bbb R$ and compose it with the function $g:\Bbb R^3\to\Bbb R^3$ that changes to spherical coordinates: $$g(r,\theta,\varphi) = (r\cos\theta\sin\varphi,r\sin\theta\sin\varphi,r\cos\varphi)$$
The result is $\tilde v(r,\theta,\varphi)=(v\circ g)(r,\theta,\varphi)$ i.e. "$v$ written in spherical coordinates". An abuse of notation is usually/almost-always commited here and we write $v(r,\theta,\varphi)$ to denote what is actually the new function $\tilde v$. I will use that notation myself now. Anyways, the chain rule states that $$\begin{pmatrix}\partial_x v & \partial_y v & \partial_z v\end{pmatrix} \begin{pmatrix} \cos\theta\sin\varphi &-r\sin\theta\sin\varphi & r\cos\theta\cos\varphi \\ \sin\theta\sin\varphi & r\cos\theta\sin\varphi &r\sin\theta\cos\varphi \\\cos\varphi & 0 & -r\sin\varphi\end{pmatrix} = \begin{pmatrix}\partial_r v & \partial_\theta v & \partial_\varphi v\end{pmatrix}$$
From this we get, for example (by inverting the matrix) that $$\partial_x = \cos\theta\sin\varphi\partial_r - \frac{\sin\theta}{r\sin\varphi}\partial_\theta + \frac{\cos\theta\cos\varphi}{r}\partial_\varphi$$
The rest will have similar expressions. Now that we know how to take partial derivatives of a real valued function whose argument is in spherical coords., we need to find out how to rewrite the value of a vector valued function in spherical coordinates. To be precise, the new basis vectors (which vary from point to point now) of $\Bbb R^3$ are found by differentiating the spherical parametrization w.r.t. its arguments (and normalizing). Thus (one example), $$\mathbf e_r = \frac{\partial_r g}{\|\partial_r g\|} = \frac{\begin{pmatrix} \cos\theta\sin\varphi & \sin\theta\sin\varphi & \cos\theta\end{pmatrix}}{\|\begin{pmatrix} \cos\theta\sin\varphi & \sin\theta\sin\varphi & \cos\theta\end{pmatrix}\|} = \begin{pmatrix} \cos\theta\sin\varphi & \sin\theta\sin\varphi & \cos\theta\end{pmatrix} \\[4ex] = \cos\theta\sin\varphi \mathbf i + \sin\theta\sin\varphi \mathbf j + \cos\theta\mathbf k$$
I don't know how to justify this without speaking of tangent spaces, but I'm sure you can ask your teacher for an explanation. After calculating the new unit vectors, you'll again have to invert the relation to obtain $\mathbf i,\mathbf j,\mathbf k$ in terms of $\mathbf e_r,\mathbf e_\theta,\mathbf e_\varphi$. But that part is just linear algebra!
Now that everything is set up, we can calculate the divergence. But what is the divergence? What I mean is, how do we write it as an abstract object that acts on functions? Here is one possibility, in terms as familiar as possible to a calculus student (there are other definitions too): $$\mathrm{div}(\cdot) = \partial_x\left(\langle \mathbf{i},\cdot\rangle\right) + \partial_y\left(\langle \mathbf{j},\cdot\rangle\right) + \partial_z\left(\langle \mathbf{k},\cdot\rangle\right)$$
Where the symbol $\langle\cdot,\cdot\rangle$ is used for the dot product. Try to convince yourself why the above formula is so.
Now just substitue all of the expressions we just derived for the basis vectors, and the differential operators. Finally, place an arbitrary vector field $$ \mathbf E = E_r(r,\theta,\varphi)\,\mathbf e_r + E_\theta(r,\theta,\varphi)\,\mathbf e_\theta + E_\varphi(r,\theta,\varphi)\,\mathbf e_\varphi$$
in place of the "$\cdot$" in the (new) expression for $\mathrm{div}$, and expand.
A: Let $E:\mathbb{R}^3\rightarrow\mathbb{R}^3$ be a vector field 
defined in Cartesian coordinates.  Let 
$M:\mathbb{R}^3\rightarrow\mathbb{R}^3$ be the mapping from spherical
coordinates to Cartesian coordinates.  Then $S=E\circ M$ is $E$ in 
spherical coordinates.  It is very easy to find the divergence of 
$E$ in terms of the partial derivatives of $S$ since 
$$
DS = (DE)(DM) \Rightarrow DE = (DS)(DM)^{-1}
$$
where $DG$ is the Jacobian matrix of vector field $G.$
Adding up the $E^1_x + E^2_y + E^3_z$ from this equation using
a computer algebra system shows that
$$
\begin{align}
\nabla\cdot E &= \big[ S^2_{\phi} \cos\phi - S^1_{\phi}\sin\phi \ + \\
& \qquad (S^3_r r + S^2_{\theta}\sin\phi + S^1_{\theta}\cos\phi)\cos\theta\sin\theta \ + \\
& \qquad \left((S^1_r\cos\phi + S^2_r\sin\phi)r - S^3_{\theta}\right)\sin^2\theta
\Big] / (r\sin\theta).
\end{align}
$$
Subscripts denote partial differentiation with 
respect to the subscript and superscripts denote vector field output
coordinate index.  The spherical coordinates convention used:
$\theta$ is the angle measured from the positive 
$z$ axis and $\phi$ is the azimuthal angle.  This formula is equivalent 
to the more common formula using $S\cdot \hat{r},$ 
$S\cdot \hat{\theta},$ and $S\cdot \hat{\phi}.$
A: https://www.cpp.edu/~ajm/materials/delsph.pdf
This PDF has all your answers and step by step approach to the Answer.
