Using Grad, Div, etc in polar coordinates

I'm trying to use the newish functions Grad and Div to work with a PDE with circular or spherical symmetry and getting what looks to me like a wrong result. However, I'm not sure of this--perhaps I just misunderstand something. Here's an example of what I want to do. I'm trying to take the gradient of a function in two dimensions that depends only on $r$.

In[42]:= Grad[V[r],{r,th},"Polar"]
Out[42]= {V'[r],0}


So, umm, this is wrong, isn't it? The gradient of a circularly symmetric function is a vector pointing towards or away from the origin. {V'[r],0} is a vector pointing along the x-axis:

In[46]:= FromPolarCoordinates[%]
Out[46]= {V'[r],0}


I can get the right answer if I do it in Cartesian coordinates and convert:

In[62]:= Grad[V[Sqrt[x^2+y^2]],{x,y}]/.(x^2+y^2)^p_.:>r^(2p)
Out[62]= {(x V'[r])/r,(y V'[r])/r}
In[63]:= ToPolarCoordinates[%]//Simplify
Out[63]= {Sqrt[((x^2+y^2) (V'[r]^2)/r^2],ArcTan[(x V'[r])/r,(y V'[r])/r]}


which, with some obvious simplifications, is {V'[r], th}, as it ought to be.

The old VectorAnalysis package handled this correctly -- or perhaps I should say, it handled it in the way I expect, but I can't get it to work now, either.

If I am misunderstanding how the new vector differentiation functions are supposed to work, I would be grateful for an explanation. Many thanks!

• The result Out[42] is a vector pointing only in r-direction, since it only has a coordinate in the first slot, i.e., only a r-coordinate (NOT an x-coordinate). The result is correct, as I see it. – Mauricio Fernández Oct 26 '15 at 14:39
• I don't think I understand you. The second coordinate of Out[42], which should be theta, the angle, is 0. Any {r, theta} with theta = 0 translates to a vector pointing along the x-axis in Cartesian coordinates. And that is in fact the result that FromPolarCoordinates gives. – Leon Avery Oct 26 '15 at 14:45
• I propose to close this question as the one based on a simple OP's mistake. – Alexei Boulbitch Oct 26 '15 at 14:59
• @LeonAvery The answer is correct as it is in e_r, e_th basis as pointed out in an answer below. The e_th component of the gradient is $1/r \partial/\partial\theta$ which vanishes for a potential that is symmetric wrt $\theta$. This is why the result by Mathematica is correct! – Lukas Oct 26 '15 at 17:56

The explanation is that the result is not quite in polar coordinates. The first and second elements of the result vector are not the length and angle of the vector. This is not how they must be interpreted, with an $r$ and $\theta$ coordinate:

Instead they are the coefficients of a radial and tangential unit vectors at the given point: $\mathbf{e}_r$ and $\mathbf{e}_\theta$, as shown in this figure:

So the result is correct, but it is not in polar coordinates. The documentation states (under Details) that

If f is a scalar, Grad[f,{x1, x2, ..., xn}, chart] returns a vector in the orthonormal basis associated with chart

• Ah! that begins to make some sense. Thank you! So, is there a straightforward way to translate this gradient vector expressed in terms of local tangents into one in the global coordinate system. (I mean, I know in this specific case it is simply a rotation by theta, but I would like a general solution such that Mathematica can do it for me for other coordinate systems as well.) – Leon Avery Oct 26 '15 at 15:19
• @LeonAvery I'm sorry, I wasn't able to find a quick, automated builtin way. You don't need to accept until you get that. – Szabolcs Oct 27 '15 at 11:03
• Thanks. I added a second answer with a non-quick, non-automated way to do it. – Leon Avery Oct 27 '15 at 11:11

The command

Grad[V[r],{r,th},"Polar"]


will give you an answer in polar coordinates, not in cartesian coordinates. That means that the results

{V'[r],0}


represents the gradient in polar coordinates

$${\rm grad} (V(r)) = V'(r) e_r + 0 e_\theta$$

If you want, you can then transform this back to cartesian coordinates, but that's usually not the point in using an alternative coordinate system, you want to describe it in the new system.

EDIT: may be this needs some background clarification. From vector/tensor analysis, depending on the transformation rule of the coordinates you have to consider a lot. In Mathematica, if you use a general function depending on all variables, you will get the general coordinates of the gradient in the respective system. Example: for polar coordinates

Grad[f[r, th], {r, th}, "Polar"]


you will get the general form for polar coordinates

{D[f[r,th],r],D[f[r,th],th]/r}


As commented in the question by Lukas, for a function NOT depending on $\theta$, the second component vanishes. This has nothing to do with the $\theta$ position of the point considered in the gradient. The gradient is a new vector, which can be represented in respect to any chosen basis. If you want to derive the gradient of a general SCALAR function $f$ $${\rm grad}f = \frac{\partial f}{\partial r} e_r + \frac{\partial f}{\partial \theta} \frac{1}{r} e_\theta$$ consider this:

Transformation of coordinates, with Cartesian coordinates $x^i$ in respect to a position INdependent orthonormal basis $\{e_i\}, e_i\cdot e_j = \delta_{ij}$ ($\delta_{ij}$ is the Kronecker symbol) and polar coordinates $\{r,\theta\} = \{u^1,u^2\}$

$$\begin{bmatrix} x^1 \\ x^2 \end{bmatrix} = u^1 \begin{bmatrix} \cos(u^2) \\ \sin(u^2) \end{bmatrix} \ , \qquad \begin{bmatrix} u^1 \\ u^2 \end{bmatrix} = \begin{bmatrix} \sqrt{(x^1)^2+(x^2)^2} \\ \arctan(x^2/x^1) \end{bmatrix}$$ The gradient of a general function $f$ is originally defined in respect to Cartesian coordinates as $${\rm grad}f = \sum_{i=1}^2\frac{\partial f}{\partial x^i}e_i$$ With the chain rule you can get this $${\rm grad}f = \sum_{i=1}^2\frac{\partial f}{\partial x^i}e_i$$= \sum_{i=1}^2\sum_{j=1}^2\frac{\partial f}{\partial u^j}\frac{\partial u^j}{\partial x^i}e_i = \sum_{j=1}^2\frac{\partial f}{\partial u^j} \sum_{i=1}^2 \frac{\partial u^j}{\partial x^i}e_i = \sum_{j=1}^2\frac{\partial f}{\partial u^j} g^j \ , at what we defined the so called gradient basis vectos $\{g^j\}$ $$g^j = \sum_{i=1}^2 \frac{\partial u^j}{\partial x^i}e_i \ .$$ The gradient basis vectors are position dependent. They are the dual vectors of the natural tangential basis vectors $g_j = \sum_{i=1}^2 \partial x^i/\partial u^j e_i$, i.e., $g_m \cdot g^n = \delta_{mn}$. For polar coordinates you get here $$g^1 = \sum_{i=1}^2 \frac{\partial u^1}{\partial x^i}e_i = \frac{x^1}{\sqrt{(x^1)^2+(x^2)^2}}e_1 + \frac{x^2}{\sqrt{(x^1)^2+(x^2)^2}}e_2 = \cos(u^2) e_1 + \sin(u^2) e_2 \ , \qquad g^2 = -\frac{\sin(u^2)}{u^1}e_1 + \frac{\cos(u^2)}{u^1}e_2 \ .$$ Sadly, these guys are, in general, nor orthogonal to each other, nor normalized. For the special case of polar coordinates, the are orthogonal, but not normalized. So if we want to describe our gradient, a vector, in respect to a orthonormal basis, we have to normalize at least the chosen basis vectors $${\rm grad}f = \sum_{j=1}^2\frac{\partial f}{\partial u^j} g^j = \sum_{j=1}^2\frac{\partial f}{\partial u^j} ||g^j|| \frac{g^j}{||g^j||} = \sum_{j=1}^2\frac{\partial f}{\partial u^j} ||g^j|| g^{*j}$$ with the normalized basis vectors $$g^{*j} = \frac{g^j}{||g^j||} \ .$$ At the end you get, that for a general function in polar coordinates, its gradient is equivalently described as $${\rm grad}f = \sum_{i=1}^2\frac{\partial f}{\partial x^i}e_i$$= \sum_{j=1}^2\frac{\partial f}{\partial u^j} ||g^j|| g^{*j} = \frac{\partial f}{\partial u^1} g^{*1} + \frac{\partial f}{\partial u^2} \frac{1}{u^1} g^{*2} = \frac{\partial f}{\partial r} e_r + \frac{\partial f}{\partial \theta} \frac{1}{r} e_\theta If you consider a function not depending on theta, $f(r,\theta) = V(r)$, then you get the result given above and what Mathematica also gives you. Mathematica does all this stuff internally, also for the gradients of higher-order tensors. You can use the gradient function of Mathematica for second-order tensors and you will get the answer using this and actually more stuff (Christoffel symbols of the second kind, ...).

This is a follow-up to @Szabolcs excellent answer, showing how to get what I (and apparently some other commenters) expected, namely an answer for the gradient in polar coordinates. The key realization is that the result returned by Grad[V[r], {r, th}, "Polar"] is in a local orthonormal basis; that is, a Cartesian coordinate system, albeit rotated. To get the gradient in polar coordinates, you need to rotate it to line up with the global Cartesian coordinate axes, then convert from Cartesian to Polar. Here's a function that does that job:

polarGrad[f_]:=
assume = CoordinateChartData["Polar", "CoordinateRangeAssumptions"][{r, th}];
rot = CoordinateTransformData["Polar"->"Cartesian", "OrthonormalBasisRotation", {r, th}];
"Cartesian"->"Polar",
];

I suspect there is a simpler way to do this with TransformedField`, but I haven't been able to figure that out yet.