It is quite possible to establish a calculative formalism (which I learned was called "dyadic formalism"), in which you may derive the correct relationships for arbitrary coordinate systems, as long as you let "differentiation" have precedence to other operations.
Let us for example, see how the Laplacian operator for 2-D polar coordinates can be derived:
We have:
\nabla=\vec{i}_{r} \frac{\partial}{\partial{r}}+\vec{i}_{\theta} \frac{\partial}{r\partial\theta}
The Laplacian operator can now be written as:
\nabla^{2}=\nabla\cdot\nabla=\vec{i}_{r}\cdot\frac{\partial}{\partial{r}}(\vec{i}_{r}\frac{\partial}{\partial{r}} + \vec{i}_{\theta}\frac{\partial}{r\partial\theta}) + \vec{i}_{\theta}\cdot\frac{\partial}{r\partial\theta}(\vec{i}_{r}\frac{\partial}{\partial{r}} + \vec{i}_{\theta}\frac{\partial}{r\partial\theta})
Now, computing all differentiations, including most importantly, those of unit vectors, we gain the correct result for the Laplacian operatior:
\nabla^{2}=\frac{\partial^{2}}{\partial{r}^{2}}+\frac{1}{r} \frac{\partial}{\partial{r}}+\frac{\partial^{2}}{r^{2}\partial\theta^{2}}
Suitably formulated, this can be used for other coordinate systems as wll, for example 3-D spherical coordinates.
When calculating the divergence or cross product of a vector, proceed in a perfectly similar way, but DO remember to differentiate those pesky unit vectors prior to the product operations; otherwise, you'll get wrong results.
----------------------------------------
Obviously, this calculation method doesn't constitute any "proof", or strict derivation on its own, but is a clever formalism that gives the correct results that, for example, the method vanhees71 describes, yields in a more rigorous manner.