Given a vector $\boldsymbol v$ on a manifold, there are an infinite number of ways to represent it as a sum of basis elements. The standard choice in differential geometry is to use the holonomic basis, that is, given the coordinate vector $\boldsymbol x=(x^1,\dots ,x^n)$, and some coordinates $q^1,\dots ,q^n$
$$\boldsymbol e_i(\boldsymbol p)= \frac{\partial\boldsymbol x}{\partial q^i}(\boldsymbol p)\tag{1}$$
And not the normalized holonomic basis,
$$\hat{\boldsymbol e}_i(\boldsymbol p)=\left\Vert\frac{\partial\boldsymbol x}{\partial q^i}(\boldsymbol p)\right\Vert^{-1}~\frac{\partial\boldsymbol x}{\partial q^i}(\boldsymbol p)\tag{2}$$
For instance, if we are using the standard basis $q^i=x^i$ then $\boldsymbol e_1=(1,\dots,0)$ etc. and in this case all of the unit vectors have a magnitude of unity, so there is no distinction between the normalized and holonomic bases.
But what happens if our coordinate vector is a more complicated function of its coordinates? E.g for 2D polar coordinates,
$$\boldsymbol x=\begin{pmatrix}x^1 \\ x^2\end{pmatrix}=\begin{pmatrix}q^1\cos(q^2) \\ q^1\sin(q^2)\end{pmatrix}$$
Then there is a distinction, and $\boldsymbol e_i\neq \hat{\boldsymbol e}_i$ in general. Therefore, it makes a difference whether we represent our vector as
$$\boldsymbol v=v^i\boldsymbol e_i\tag{3}$$
Or as
$$\boldsymbol v=\hat{v}^i\hat{\boldsymbol e}_i\tag{4}$$
Using $(1),(2)$, the different bases are related by
$$\hat{\boldsymbol e}_i=\frac{1}{\Vert\boldsymbol e_i\Vert}\boldsymbol e_i$$
Which, using $(3),(4)$ means the normalized and non-normalized components are related by
$$ v^i=\frac{1}{\Vert\boldsymbol e_i\Vert}\hat{v}^i\tag{*}$$
So, given a vector $\boldsymbol v$ with its coordinates $v^i$ represented in the holonomic basis, the formula for its covariant derivative is
$$(\nabla v)^i{}_j=\partial_jv^i+\Gamma^i_{jk}v^k$$
The fact that the formula Wolfram Mathworld gives has that weird factor of $\frac{1}{g_{ii}}$ in front should already tip you off that something funny is going on. Doing some algebra, this gives the $(1,2)$ or $(r,\theta)$ component of the covariant derivative of $\boldsymbol v$ as
$$(\nabla v)^r{}_\theta=\partial_{\theta}v^r-rv^\theta$$
However, if we want to represent our results with the normalized components $\hat{v}^i$, this becomes
$$(\nabla v)^r{}_\theta=\partial_{\theta}\left(\frac{1}{\Vert\boldsymbol e_r\Vert}\hat{v}^r\right)-r~\frac{1}{\Vert\boldsymbol e_\theta\Vert}\hat{v}^\theta$$
Given that the coordinate vector is $\boldsymbol x=(q^1\cos(q^2),q^1\sin(q^2),q^3)=(r\cos\theta,r\sin\theta,z)$ you can calculate by yourself that
$$\Vert\boldsymbol e_r\Vert=1~~,~~\Vert\boldsymbol e_\theta\Vert=r$$
Using that result, and $(*)$, we get
$$(\nabla v)^r{}_{\theta}=\partial_\theta \hat{v}^r-\hat{v}^\theta$$
But, hm.... something still isn't right. There appears to be a missing factor of $1/r$. Then we realize - the tensor components on the left hand side are still being represented with the holonomic basis! Recall that we can represent a $(1,1)$ tensor as
$$\mathbf T=T^i{}_j~\boldsymbol e_i\otimes \boldsymbol e^j$$
Where $\boldsymbol e^j=g^{jk}\boldsymbol e_k$ is the dual basis. Once again we can make a dual normalized basis,
$$\hat{\boldsymbol e}^j=\frac{1}{\Vert\boldsymbol e^j\Vert}\boldsymbol e^j$$
So if we represent the tensor $\mathbf T$ with its normalized components,
$$\mathbf T=\hat{T}^i{}_j~\hat{\boldsymbol e}_i\otimes \hat{\boldsymbol e}^j$$
They are related to the holonomic components $T^i{}_j$ via
$$T^i{}_j=\frac{1}{\Vert\boldsymbol e_i\Vert}\frac{1}{\Vert\boldsymbol e^j\Vert}\hat{T}^i{}_j\tag{**}$$
You can compute the dual basis vectors for cylindrical coordinates to be
$$\boldsymbol e^r=\begin{pmatrix}\cos \theta & \sin\theta & 0\end{pmatrix} \\
\boldsymbol e^\theta=\frac{1}{r^2}\begin{pmatrix}-r\sin\theta & r\cos\theta & 0\end{pmatrix} \\ \boldsymbol e^z=\begin{pmatrix}0&0&1\end{pmatrix}$$
Hence
$$\Vert\boldsymbol e^\theta\Vert=\frac{1}{r}$$
So using that result, and $(**)$, you get
$$(\nabla v)^{r}{}_\theta=\frac{1}{\Vert\boldsymbol e_r\Vert}\frac{1}{\Vert\boldsymbol e^\theta\Vert}(\widehat{\nabla v})^r{}_\theta=r~(\widehat{\nabla v})^r{}_\theta$$
And hence
$$\boxed{(\widehat{\nabla v})^r{}_\theta=\frac{1}{r}(\nabla v)^{r}{}_\theta=\frac{1}{r}\partial_\theta v^r-\frac{1}{r}v^\theta}$$
The formula that you were looking for.
MORAL OF THE STORY:
All of the differential geometry formulas in almost any reputable textbook are given using the holonomic basis. (Why? Because they're generally simpler.) So - stick to the holonomic basis for as long as you possibly can. If you must convert to the orthonormal basis, do it at the end, after you have finished all of your calculations. If you do not, you will use formulas that work when the components are expressed in a holonomic basis, on a vector where you have represented the components using a normalized basis, and you are at a high risk of getting incorrect results.
Reference - Christoffel symbols and covariant derivative components in cylindrical coordinates
Christoffel symbols
$$\Gamma^r_{\theta\theta}=-r \\ \Gamma^{\theta}_{r\theta}=\Gamma^\theta_{\theta r}=\frac{1}{r}\\\text{all others zero}$$
Covariant derivative
$$\nabla\boldsymbol u=\begin{pmatrix}\partial_r u^r & \partial_\theta u^r-ru^\theta & \partial_z u^r \\ \partial_r u^\theta+\frac{1}{r}u^\theta & \partial_\theta u^\theta +\frac{1}{r}u^r & \partial_z u^\theta \\ \partial_r u^z & \partial_\theta u^z & \partial_z u^z\end{pmatrix} \\ \text{(holonomic basis)}$$
$$\nabla \boldsymbol u=\begin{bmatrix}
\partial_ru^r & \frac{1}{r}\left(\partial_\theta u^r-u^\theta\right) & \partial_z u^r \\
\partial_r u^\theta & \frac{1}{r}\left(\partial_\theta u^\theta+u^r\right) & \partial_z u^\theta \\
\partial_r u^z & \frac{1}{r}\partial_\theta u^z & \partial_z u^z
\end{bmatrix} \\ \text{(normalized basis)}$$
The second matrix agrees with all of the formulas given on Wolfram Mathworld. It can be obtained by substituting $u^\theta\to u^\theta/r$ in the expression given in the holonomic basis and then multiplying through by the appropriate magnitudes of the basis elements.
Given below is a Mathematica script that will generate both of these matrices.
$Assumptions = (_ \[Element] Reals &&
r > 0); (*Tells Mathematica to assume all variables are real, and that the radius is positive*)
dim = 3;
q = {r, \[Theta], z};
p = {r Cos[\[Theta]], r Sin[\[Theta]], z};
e = Table[D[p, q[[i]]], {i, 1, dim}];
norme = Simplify[Table[Norm[e[[i]]], {i, 1, dim}]];
g = Simplify[
Table[D[p, q[[i]]] . D[p, q[[j]]], {i, 1, dim}, {j, 1, dim}]];
inverseg = Simplify[Inverse[g]];
duale = Table[Sum[inverseg[[i, j]]*e[[j]], {i, 1, dim}], {j, 1, dim}];
normduale = Simplify[Table[Norm[duale[[j]]], {j, 1, dim}]];
\[CapitalGamma] =
Table[Sum[
inverseg[[k, m]]/
2 (D[g[[i, k]], q[[j]]] + D[g[[j, k]], q[[i]]] -
D[g[[i, j]], q[[k]]]), {k, 1, dim}], {m, 1, dim}, {i, 1,
dim}, {j, 1, dim}];
u = {ur[r, \[Theta], z], u\[Theta][r, \[Theta], z],
uz[r, \[Theta], z]};
uNORMALIZED = Table[1/norme[[i]]*u[[i]], {i, 1, dim}];
Delu = Table[(D[u[[i]], q[[j]]]) +
Sum[\[CapitalGamma][[i, j, k]] u[[k]], {k, 1, dim}], {i, 1,
dim}, {j, 1, dim}]
DeluNORMALIZED =
Table[norme[[i]]*
normduale[[j]]*((D[uNORMALIZED[[i]], q[[j]]]) +
Sum[\[CapitalGamma][[i, j, k]] uNORMALIZED[[k]], {k, 1,
dim}]), {i, 1, dim}, {j, 1, dim}]
And here it is in action:
