26
$\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]
Improve this question
$\endgroup$
0

4 Answers 4

Reset to default
23
$\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.

Improve this answer
$\endgroup$
7
  • $\begingroup$ (+1) Thank you very much jose. :) Can you also elaborate on what HodgeDual does? :) $\endgroup$ Commented Feb 13, 2019 at 9:46
  • 4
    $\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$ Commented Feb 13, 2019 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$ Commented Feb 13, 2019 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$ Commented Feb 13, 2019 at 10:31
  • 7
    $\begingroup$ +1 Crystal clear! For readers seeking more math, here is a related question on math.SE. $\endgroup$ Commented Feb 14, 2019 at 17:29
11
$\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 *)
Improve this answer
$\endgroup$
2
  • $\begingroup$ Comments are not for extended discussion; this conversation has been moved to chat. $\endgroup$ Commented Feb 13, 2019 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$ Commented Feb 13, 2019 at 11:03
5
$\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 Properties & 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.

Improve this answer
$\endgroup$
0
4
$\begingroup$

I am limiting myself to fields defined in the three dimensional Euclidean point space. The curl of a tensor can be found in these simple steps:

  1. Take the simple composition of the second-order tensor, T, with the LeviCivitaTensor[3]. This is effected by the command, 

    Dot[T, LeviCivitaTensor[3]]

  2. Take the transpose of the divergence of the above result using the command:

    curl[T_] := Transpose[Div[Dot[T, LeviCivitaTensor[3]], {x, y, z}]]

The result you get here is the second-order tensor that is the curl of the tensor T.

Improve this answer
$\endgroup$
8
  • $\begingroup$ (+1) Hi there, and welcome to this community. :) Have you tested your formula against a tensor of an arbitrary rank n and in different curivilinear coordinates such as cylindrical or spherical? Your answer is somehow similar to that of xzczd. $\endgroup$ Commented May 17, 2019 at 14:36
  • $\begingroup$ Thanks for your question. $\endgroup$ Commented May 17, 2019 at 17:44
  • $\begingroup$ To H.R. 1. Formula tested correctly in Cylindrical and Spherical Polar coordinates using the same example problem in Mathematica 12, transformed to cylindrical and spherical polar coordinates; 2. Only interested in 3D Euclidean Point Space. 3. No, this method is different from xzczd. Here, we use transpose of the divergence of the contraction of the tensor argument with the LeviCivitaTensor[3]. I have proved elsewhere that it should give the same result but from different angles except for the scalar multiplier of 1/2. $\endgroup$ Commented May 17, 2019 at 20:04
  • $\begingroup$ @H.R. If I understand it correctly, Fakinlede's answer is just for the 2nd part of the question i.e. "How it can be related to the definition given in $(2)$?" Then I think the first step should be modified to Dot[Transpose@T, LeviCivitaTensor[3]]. ( @OAFakinlede Just compare your result to the CurlTen in the question. ) $\endgroup$ Commented May 18, 2019 at 3:54
  • $\begingroup$ OK, that's fine. But one point I would like to mention is that you can have a tensor of arbitrary rank n in a 3D Euclidean space. Indeed, we encounter a lot of such creatures in mathematical physics. :) $\endgroup$ Commented May 18, 2019 at 7:08

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.