20
$\begingroup$

I have a usual mathematical background in vector and tensor calculus. I was trying to use the differential operators of Mathematica, namely Grad, Div and Curl. According to my knowledge, the definitions of Mathematica for Grad and Div coincides with those usually employed in tensor calculus, that is to say

\begin{align*} \text{grad}\mathbf{T}&:=\sum_{k=1}^{3}\frac{\partial\mathbf{T}}{\partial x_k}\otimes \mathbf{e}_k\\ \text{div}\mathbf{T}&:=\sum_{k=1}^{3}\frac{\partial\mathbf{T}}{\partial x_k}\cdot\mathbf{e}_k \\ \tag{1} \end{align*}

for any tensor $\mathbf{T}$ of rank $n\ge1$. $x_k$'s are Cartesian coordinates and $\mathbf{e}_i$'s are the standard basis for $\mathbb{R}^3$. $\otimes$ and $\cdot$ are the usual generalized outer and inner products which are also defined in Mathematica by Outer and Inner. The usual definition that I know from tensor calculus for the Curl is as follows \begin{align*} \text{curl}\mathbf{T}&:=\sum_{k=1}^{3}\mathbf{e}_k\times\frac{\partial\mathbf{T}}{\partial x_k}. \tag{2} \end{align*} However, it turns out that Mathematica's definition for curl is totally different. For example, it returns the Curl of a second order tensor as a scalar, while according to $(2)$ it should be a second order tensor.

I couldn't find a precise definition of Mathematica for Curl in the documents. I am wondering what this definition is. What is the motivation for this? and How it can be related to the definition given in $(2)$?

Below is a simple piece of code for you to observe the outputs of Mathematica when we apply the Grad, Div and Curl operators to scalar, vector and second order tensor fields. I would like to draw your attention to some observations. Curl of a scalar is returned as a second order tensor, which I don't understand why! Curl of a vector coincides with the usual definition of Curl used in vector calculus. Curl of second order tensor is returned as a scalar, which I don't understand again.

Var={Subscript[x, 1],Subscript[x, 2],Subscript[x, 3]};
Sca=\[Phi][Subscript[x, 1],Subscript[x, 2],Subscript[x, 3]];
Vec={Subscript[v, 1][Subscript[x, 1],Subscript[x, 2],Subscript[x, 3]],Subscript[v, 2][Subscript[x, 1],Subscript[x, 2],Subscript[x, 3]],Subscript[v, 3][Subscript[x, 1],Subscript[x, 2],Subscript[x, 3]]};
Ten=Table[Subscript[T, i,j][Subscript[x, 1],Subscript[x, 2],Subscript[x, 3]],{i,1,3},{j,1,3}];
MatrixForm[Grad[Sca, Var]]
MatrixForm[Grad[Vec, Var]]
MatrixForm[Grad[Ten, Var]]
MatrixForm[Div[Sca, Var]]
MatrixForm[Div[Vec, Var]]
MatrixForm[Div[Ten, Var]]
MatrixForm[Curl[Sca, Var]]
MatrixForm[Curl[Vec, Var]]
MatrixForm[Curl[Ten, Var]]

I will be happy if someone can reproduce the following result for the curl of a second order tensor with Mathematica's Curl function.

\begin{align*} \text{curl}\mathbf{T}&:=\sum_{k=1}^{3}\mathbf{e}_k\times\frac{\partial\mathbf{T}}{\partial x_k}=\sum_{k=1}^{3}\mathbf{e}_k\times\frac{\partial}{\partial x_k}\left(\sum_{i=1}^{3}\sum_{j=1}^{3}T_{ij}\mathbf{e}_i\otimes\mathbf{e}_j\right)\\ &=\sum_{k=1}^{3}\sum_{i=1}^{3}\sum_{j=1}^{3}\frac{\partial T_{ij}}{\partial x_k}(\mathbf{e}_k\times\mathbf{e}_i)\otimes\mathbf{e}_j\\ &=\sum_{k=1}^{3}\sum_{i=1}^{3}\sum_{j=1}^{3}\sum_{m=1}^{3}\epsilon_{kim}\frac{\partial T_{ij}}{\partial x_k}\mathbf{e}_m\otimes\mathbf{e}_j \tag{3} \end{align*}

where $\epsilon_{kim}$ is the LeviCivitaTensor for $3$ dimensions. Consequently, we get

\begin{align*} \left(\text{curl}\mathbf{T}\right)_{mj}=\sum_{k=1}^{3}\sum_{i=1}^{3}\epsilon_{kim}\frac{\partial T_{ij}}{\partial x_k}. \tag{4} \end{align*}

Implementing $(4)$ in Mathematica, we obtain

CurlTen = Table[
   Sum[
    LeviCivitaTensor[3][[k, i, m]] 
     D[Subscript[T, i, j][Subscript[x, 1], Subscript[x, 2], Subscript[x, 3]], {Subscript[x, k]}], {k, 1, 3}, {i, 1, 3}],
   {m, 1, 3}, {j, 1, 3}];
MatrixForm[CurlTen]
$\endgroup$
17
$\begingroup$

The definition used (motivated by exterior calculus) is as follows:

Given a rectangular array $a$ of depth $n$, with dimensions $\{d, ..., d\}$ (so there are $n$ $d$'s) and a list $x = \{x_1, ..., x_d\}$ of variables, then

Curl[a, x] == (-1)^n (n+1) HodgeDual[Grad[a, x], d]

If $a$ has depth $n$, then Grad[a, x] has depth $n+1$, and therefore HodgeDual[Grad[a, x], d] has depth $d-(n+1)$. Clearly then we need $n < d$. Note that $n = 0$ is admitted, i.e. we can take the curl of a scalar function.

In the traditional case of $d=3$ and $n=1$ we have $d-(n+1)=1$ and that's why the curl of a vector is also a vector.

The HodgeDual operation starts by antisymmetrizing its first argument, and hence implicitly we really have something like

Curl[a, x] == (-1)^n (n+1) HodgeDual[Symmetrize[Grad[a, x], Antisymmetric[All]], d]

Finally, one more comment: The definition given assumes that we work with tensors given in components in a Cartesian coordinated and orthonormal basis. If that wasn't the case we would have to insert some additional metric factors. That's handled by the third argument of Curl for a variety of alternative coordinate systems.

$\endgroup$
  • $\begingroup$ (+1) Thank you very much jose. :) Can you also elaborate on what HodgeDual does? :) $\endgroup$ – H. R. Feb 13 at 9:46
  • 1
    $\begingroup$ HodgeDual is an algebraic operation. It does not contain derivatives. Note for example that Cross of d-1 vectors in dimension d is the HodgeDual of the tensor product of those vectors, modulo a factor (d-1)! due to antisymmetrization. Roughly speaking, HodgeDual takes some directions in space and gives you the orthogonal directions. It's just that rather than giving you lists of vectors, it gives you antisymmetric tensors, the forms of exterior algebra in vector spaces. $\endgroup$ – jose Feb 13 at 9:59
  • 1
    $\begingroup$ Well, I think that I need to read some math! :D I am not familiar with exterior algebra and cross of many vectors. Do you have any simple reference in mind? :) $\endgroup$ – H. R. Feb 13 at 10:02
  • 3
    $\begingroup$ @H.R. The HodgeDual or Hodge star operator gives you the complementary subspace that 'completes' the vector space. There is a nice introduction to it (and exterior calculus in general) in Keenan Cranes book Discrete Differential Geometry: An Applied Introduction in chapter four. $\endgroup$ – Thies Heidecke Feb 13 at 10:31
  • 7
    $\begingroup$ +1 Crystal clear! For readers seeking more math, here is a related question on math.SE. $\endgroup$ – Silvia Feb 14 at 17:29
7
$\begingroup$

I fail to find a reference for the definition used by Curl, but manage to figure out how the Curl is defined.

I'll use Einstein summation convention for simplicity. The defintion used by Curl is summarized as follows:

  1. Curl of a scalar / zero-order tensor is defined as

$$\text{curl}\ \phi= \epsilon_{kim}\frac{\partial\phi}{\partial x_m}$$

  1. Curl of a vector / first-order tensor is defined as

$$\text{curl}\ \mathbf{v}= \epsilon_{kim}\frac{\partial v_m}{\partial x_i}$$

  1. Curl of a matrix / second-order tensor is defined as

$$\text{curl}\ \mathbf{T}=\color{red}{\frac{1}{2}} \epsilon_{kim}\frac{\partial T_{im}}{\partial x_k}$$

Check

ϵ = LeviCivitaTensor[3];

Table[Sum[ϵ[[k, i, m]] D[Sca, Var[[m]]], {m, 3}], {k, 3}, {i, 3}] == 
ϵ.D[Sca, {Var}] == 
TensorContract[ϵ\[TensorProduct]D[Sca, {Var}], {{3, 4}}] == 
Curl[Sca, Var]
(*True*)

Table[Sum[ϵ[[k, i, m]] D[Vec[[m]], Var[[i]]], {m, 3}, {i, 3}], {k, 3}] == 
TensorContract[ϵ\[TensorProduct]D[Vec, {Var}], {{3, 4}, {2, 5}}] == 
Curl[Vec, Var]
(*True*)

1/2 Sum[ϵ[[k, i, m]] D[Ten[[i, m]], Var[[k]]], {k, 3}, {i, 3}, {m, 3}] == 
1/2 TensorContract[ϵ\[TensorProduct]D[Ten, {Var}], {{2, 4}, {3, 5}, {1, 6}}] == 
Curl[Ten, Var]
(*True*)

The second question is simple:

Table[Sum[ϵ[[k, i, m]] D[Ten[[i, j]], Var[[k]]], {k, 3}, {i, 3}], {m, 3}, {j, 3}] == 
TensorContract[ϵ\[TensorProduct]D[Ten, {Var}], {{2, 4}, {1, 6}}] == 
(Curl[#, Var] & /@ Transpose@Ten // Transpose)
(* True *)
$\endgroup$
  • $\begingroup$ Comments are not for extended discussion; this conversation has been moved to chat. $\endgroup$ – Kuba Feb 13 at 10:40
  • 2
    $\begingroup$ To the downvoter, I am interested in what's wrong with my question, would you please elaborate. I'm not trying to complain here, I just want improve my answer if possible. $\endgroup$ – xzczd Feb 13 at 11:03
2
$\begingroup$

I suggest to modify your simple piece of code to make the output more informative. I define a function dim[] which is similar to the Dimensions[] function but is restricted to matrices and arrays. For example, try Dimensions[x+y] which returns {2} but dim[x+y] returns 0.

ClearAll[dim, id];
dim[x : (_List | _StructuredArray)] := Dimensions@x;
dim[x_Div] = $Failed;
dim[x_] = 0;
id[x_] := (Print[dim[x // Expand]]; x // MatrixForm);

Print["Grad"];
id@Grad[Sca, Var]
id@Grad[Vec, Var]
id@Grad[Ten, Var]
Print["Div"];
id@Div[Sca, Var]
id@Div[Vec, Var]
id@Div[Ten, Var]
Print["Curl"];
id@Curl[Sca, Var]
id@Curl[Vec, Var]
id@Curl[Ten, Var]

You wrote

I couldn't find a precise definition of Mathematica for Curl in the documents

The documentation for Curl[] is long and typical of of some of the more obscure Mathematica function documentation. The answer is buried in Properites & Relations where it states:

Curl[f, {x, y, z}] == (-1)^r (r + 1) HodgeDual[ Symmetrize[Grad[f, {x, y, z}], Antisymmetric[All]]]

but HodgeDual[] is even more obscure so I am not surprised that the precise definition is not given. It would be a good idea if it was available to users. So the answer to your first question is that it is still obscure. Even with a definition, one usually tries various examples and examines the results. Sometimes you have to use the function in not obvious ways. In any case, if a function does not do what you want, you can attempt to write your own function that does exactly what you want. You may decide that this is the better option in your situation.

It is good that you seem to have implemented your own version. I suggest an alternative version and a simple test

CurlTen == Transpose@Array[Curl[Ten[[All, #]], Var] &, 3]

which returns True. I hope this helps you in some small way even though I can't answer your first question.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.