Message Boards Message Boards

Coordinate-free Tensorial/Vectorial Calculus

GROUPS:

There are a lot of Mathematica packages for vector/tensor calculus: Tensorial, Advanced Tensor Analysis, Ricci, TensoriaCalc, grt, xAct are just a few mentions.

The main thing I don't like about these packages is that the majority of them (if not all) are not up to date with the current Mathematica version. Some of them haven't seen an update in decades. The other major drawback of these packages is the cumbersome notation, declaration and the fact that almost all are to be used with General Relativity in mind and use coordinate notation. The results are not really coordinate-free.

With this in mind, I have developed a "Package" that can deal with symbolic vectors/tensors in a coordinate-free form and the package is human-readable (some packages implementations in this matter are simply indecipherable).

In the first part, I'll present the code with simple explanations and the second part some simple examples.

These functions are needed to clear the OenOwnValues, DownValues, UpValues, SubValues. Mathematica has no built-in way of doing it.

(* Prevent evaluation *)
SetAttributes[{ClearOwnValues, ClearDownValues, ClearUpValues, ClearSubValues}, HoldFirst]

(* Always return true for commodity later *)
ClearOwnValues[var_Symbol]  := (OwnValues@var  = {}; True)
ClearDownValues[var_Symbol] := (DownValues@var = {}; True)
ClearUpValues[var_Symbol]   := (UpValues@var   = {}; True)
ClearSubValues[var_Symbol]  := (SubValues@var  = {}; True)

(* Delete Values of "f" if they match the input *)
ClearDownValues[expr:f_Symbol[___]]     := (DownValues@f = DeleteCases[DownValues@f, _?(!FreeQ[First@#, HoldPattern@expr] &)]; True)
ClearUpValues[expr:f_Symbol[___]]       := (UpValues@f   = DeleteCases[UpValues@f,   _?(!FreeQ[First@#, HoldPattern@expr] &)]; True)
ClearSubValues[expr:f_Symbol[___][___]] := (SubValues@f  = DeleteCases[SubValues@f,  _?(!FreeQ[First@#, HoldPattern@expr] &)]; True)

The Define-family functions are the main core of this package, they are used to define all kinds of relationship and can be easily expanded. You can define a Symbol "var" or a function "fun[var, ___]" to be a certain "type". "var" cannot have OwnValues, otherwise it will be evaluated. You can only define symbols. the function "fun" should be primarily used as script-like functions as Subscript, SuperHat, etc.

(* Prevent evaluation *)
SetAttributes[Define$Internal, HoldAll]

(* 
    Internal version of Define, all others versions are a call of this function.
    The OwnValues and/or DownValues are cleared immediately, this step is to avoid rule-parsing.
    The possible types are: Real, Imaginary, Tensor and Constant. For Tensor type additional parameters are needed: rank and dimension.
*)
Define$Internal[(var_Symbol /; ClearOwnValues@var) | (var:Except[Hold, fun_][head_Symbol /; ClearOwnValues@head, ___] /; ClearDownValues@var),
    type:"Real" | "Imaginary" | "Tensor" | "Constant",
    rank_Integer:2, dim_Integer:3] := Module[{tag},

    (* All expressions are defined as UpValues, assign it to the corresponding tag *)
    (* UpValues cannot be deeply nested, hence the need to assign it to the "tag" *)
    tag = If[Head@var === Symbol, var, fun];

    Which[
       type === "Real", (* Typical properties needed for real quantities *)
         Evaluate@tag /: Element[var, Reals] = True;
         Evaluate@tag /: Re[v:var] := v;
         Evaluate@tag /: Im[v:var] := 0;
         Evaluate@tag /: Conjugate[v:var] := v;
         Evaluate@tag /: Abs[v:var] := RealAbs@v;
       ,        
       type === "Imaginary", (* Typical properties needed for Imaginary quantities *)
         Evaluate@tag /: Element[var, Reals] = False;
         Evaluate@tag /: Re[v:var] = 0;
         Evaluate@tag /: Im[v:var] := v/I;
         Evaluate@tag /: Conjugate[v:var] := -v;
         Evaluate@tag /: Abs[v:var] := RealAbs@Im@v;
       ,        
       type === "Tensor", (* For compativility with Mathematica current Tensor-functions *)
         Evaluate@tag /: ArrayQ@var = rank != 0;
         Evaluate@tag /: TensorQ@var = rank != 0;
         Evaluate@tag /: MatrixQ@var = rank == 2;
         Evaluate@tag /: VectorQ@var = rank == 1;
         Evaluate@tag /: ListQ@var = rank != 0;
         Evaluate@tag /: ScalarQ@var = rank == 0;
         Evaluate@tag /: TensorRank@var = rank;
         Evaluate@tag /: TensorDimensions@var = ConstantArray[dim, {rank}];
         Evaluate@tag /: Element[var, Arrays@TensorDimensions@var] = True;
         Evaluate@tag /: Element[var, Matrices@{dim, dim}] = rank == 2;
         Evaluate@tag /: Element[var, Vectors@dim] = rank == 1;
       ,        
       type === "Constant", (* A constant has zero "derivative" *)
         Evaluate@tag /: ConstantQ@var = True;
         Evaluate@tag /: grad@var = 0;
         Evaluate@tag /: div@var = 0;
         Evaluate@tag /: curl@var = 0;
         Evaluate@tag /: DotNabla[_, var] = 0;
         Evaluate@tag /: D[var, __] = 0;  
         Evaluate@tag /: Dp[var, __] = 0; 
         Evaluate@tag /: Dt[var, ___] = 0;
         Evaluate@tag /: Delta[var] = 0;
       ,        
       True, $Failed]    
]

(* Assign more than one variable *)
Define$Internal[vars__ /; Length@{vars} > 1, type:"Real" | "Imaginary" | "Tensor" | "Constant", rank_Integer:2, dim_Integer:3] :=(
    Define$Internal[#, type, rank, dim] & /@ Hold /@ Hold@vars // ReleaseHold;) (* Hacky-way of passing Hold down *)
Define$Internal[Hold@var_, type:"Real" | "Imaginary" | "Tensor" | "Constant", rank_Integer:2, dim_Integer:3] := Define$Internal[var, type, rank, dim]

(* Main Define functions *)
SetAttributes[{DefineReal, DefineImaginary, DefineTensor, DefineConstant}, HoldAll]
DefineReal[vars__] := Define$Internal[vars, "Real"]
DefineImaginary[vars__] := Define$Internal[vars, "Imaginary"]
DefineTensor[vars__, rank_Integer:2, dim_Integer:3] := Define$Internal[vars, "Tensor", rank, dim]
DefineConstant[vars__] := Define$Internal[vars, "Constant"]

(* Define multiple things at once *)
SetAttributes[{DefineRealTensor, DefineConstantTensor, DefineRealConstantTensor}, HoldAll]
DefineRealTensor[vars__, rank_Integer:2, dim_Integer:3] := (DefineReal@vars; DefineTensor[vars, rank, dim];)
DefineConstantTensor[vars__, rank_Integer:2, dim_Integer:3] := (DefineConstant@vars; DefineTensor[vars, rank, dim];)
DefineRealConstantTensor[vars__, rank_Integer:2, dim_Integer:3] := (DefineReal@vars; DefineConstant@vars; DefineTensor[vars, rank, dim];)

Now it is possible to define tensorial variables and make them behave as tensor with current Mathematica implementation. Some built-in functions needed to be redefined to work with symbolic tensors. An example of this necessity is:

(* Define two tensors a and b *)
DefineTensor[a, b, 2]
TensorRank[2*a - 3*b] (* Return 2 *)
TensorQ[a] (* Return True *)
TensorQ[2*a] (* Return False *)

Mathematica function TensorQ don't know that a scalar times a tensor is a tensor. The following code is for refifinition:

Unprotect[TensorQ, VectorQ, TensorRank, Dot, Cross, TensorProduct]

(* Numbers are always scalar/constant. These functions are not built-in. *)
ScalarQ[a_?NumericQ]   := True
ConstantQ[a_?NumericQ] := True

(* Complexes *)
   TensorQ[(Re|Im|Conjugate)[a_]] := TensorQ@a
   VectorQ[(Re|Im|Conjugate)[a_]] := VectorQ@a
   ScalarQ[(Re|Im|Conjugate)[a_]] := ScalarQ@a
 ConstantQ[(Re|Im|Conjugate)[a_]] := ConstantQ@a
TensorRank[(Re|Im|Conjugate)[a_]] := TensorRank@a

(* Plus *)
TensorQ[(a_?TensorQ) + (b_?TensorQ)] := TensorRank@a === TensorRank@b
VectorQ[(a_?VectorQ) + (b_?VectorQ)] := True
ScalarQ[(a_?ScalarQ) + (b_?ScalarQ)] := True
ConstantQ[(a_?ConstantQ) + (b_?ConstantQ)] := True

(* Times *)
TensorQ[(a__?ScalarQ) * (b_?TensorQ)] := True
VectorQ[(a__?ScalarQ) * (b_?VectorQ)] := True
ScalarQ[(a__?ScalarQ) * (b_?ScalarQ)] := True
ConstantQ[(a_?ConstantQ /; ScalarQ@a) * (b_?ConstantQ)] := True
(* Pass scalars out of Dot and Cross, as is done in TensorProduct *)
Dot[a___, Times[b_, s__?ScalarQ], c___] := Times[s, Dot[a, b, c]]
Cross[a_, Times[b_, s__?ScalarQ]] := Times[s, Cross[a, b]]
Cross[Times[a_, s__?ScalarQ], b_] := Times[s, Cross[a, b]]

(* Dot *)
TensorQ[(a_?TensorQ) . (b_?TensorQ)] /; TensorRank@a + TensorRank@b - 2 >= 1 := True
VectorQ[(a_?TensorQ) . (b_?TensorQ)] /; TensorRank@a + TensorRank@b - 2 == 1 := True
ScalarQ[(a_?VectorQ) . (b_?VectorQ)] := True
ConstantQ[(a_?ConstantQ /; TensorQ@a) . (b_?ConstantQ /; TensorQ@b)] := True
(* Automatically evaluate to zero, as TensorProduct *)
Dot[a___, 0, b___] := 0

(* Cross *)
TensorQ[(a_?VectorQ) \[Cross] (b_?VectorQ)] := True
VectorQ[(a_?VectorQ) \[Cross] (b_?VectorQ)] := True
ConstantQ[(a_?ConstantQ /; VectorQ@a) \[Cross] (b_?ConstantQ /; VectorQ@b)] := True
(* Perpendicular vectors automatically evalute to zero *)
Cross[a_?VectorQ, a_?VectorQ] := 0
(* Automatically evaluate to zero, as TensorProduct *)
Cross[a___, 0, b___] := 0
(* Return single argument as Dot, Times and TensorProduct *)
Cross[a_] := a

(* Tensor Product *)
TensorQ[(a_?TensorQ) \[TensorProduct] (b_?TensorQ)] := True
ConstantQ[(a_?ConstantQ /; TensorQ@a) \[TensorProduct] (b_?ConstantQ /; TensorQ@b)] := True

(* Power *)
ScalarQ@Power[a_?ScalarQ, b_?ScalarQ] := True
ScalarQ[1/a_?ScalarQ] := True
ConstantQ@Power[a_?ConstantQ /; ScalarQ@a, b_?ConstantQ /; ScalarQ@b] := True
ConstantQ[1/a_?ConstantQ /; ScalarQ@a] := True

(* grad *)
grad[_?ConstantQ] := 0
TensorQ@grad[a_?ScalarQ] := True
VectorQ@grad[a_?ScalarQ] := True
TensorQ@grad[a_?TensorQ] := True
TensorRank@grad[a_?ScalarQ] := 1
TensorRank@grad[a_?TensorQ] := TensorRank@a + 1

(* div *)
div[_?ConstantQ] := 0
TensorQ@div[a_?TensorQ /; TensorRank@a >= 2] := True
VectorQ@div[a_?TensorQ /; TensorRank@a == 2] := True
ScalarQ@div[a_?VectorQ] := True
TensorRank@div[a_?TensorQ] := TensorRank@a - 1

(* curl *)
curl[_?ConstantQ] := 0
TensorQ@curl[a_?VectorQ] := True
VectorQ@curl[a_?VectorQ] := True
TensorRank@curl[a_?VectorQ] := 1

(* DotNabla *)
DotNabla[_, _?ConstantQ] := 0

(* Dp *)
Dp[_?ConstantQ, args__] := 0
TensorQ@Dp[a_, args__] := TensorQ@a
VectorQ@Dp[a_, args__] := VectorQ@a
ScalarQ@Dp[a_, args__] := ScalarQ@a
TensorRank@Dp[a_, args__] := TensorRank@a

(* Delta *)
Delta[_?ConstantQ] := 0
TensorQ@Delta[a_] := TensorQ@a
VectorQ@Delta[a_] := VectorQ@a
ScalarQ@Delta[a_] := ScalarQ@a
TensorRank@Delta[a_] := TensorRank@a

(* List *)
Dp[a_List, args__] := Dp[#, args] & /@ a

(* Don't assume anything is a scalar/constant *)
ScalarQ[a_]   := False
ConstantQ[a_] := False

Protect[TensorQ, VectorQ, TensorRank, Dot, Cross, TensorProduct]

Where we have defined the Tensor-functions: grad, div, curl which are self-explanatory; Dp is the partial derivative, Delta gives the variation of a quantity, somewhat related to Dp, and DotNabla is (for the lack of better name) the convective derivative.

For better print, we'll define the following notation:

(* Hacky-way to create parenthesis *)
MakeBoxes[Parenthesis[a_], _] := MakeBoxes[a.1][[1, 1]]
MakeBoxes[grad[a_], form:TraditionalForm] := TemplateBox[{MakeBoxes@Parenthesis@a}, "grad", Tooltip -> Automatic,
    DisplayFunction :> (RowBox@{"\[Del]", #1} &)]
MakeBoxes[div[a_], form:TraditionalForm] := TemplateBox[{MakeBoxes@Parenthesis@a}, "div", Tooltip -> Automatic,
    DisplayFunction :> (RowBox@{"\[Del]\[CenterDot]", #1} &)]
MakeBoxes[curl[a_], form:TraditionalForm] := TemplateBox[{MakeBoxes@Parenthesis@a}, "curl", Tooltip -> Automatic,
    DisplayFunction :> (RowBox@{"\[Del]\[Cross]", #1} &)]
MakeBoxes[DotNabla[a_, b_], form:TraditionalForm] := TemplateBox[{MakeBoxes@Parenthesis@a, MakeBoxes@Parenthesis@b}, "DotNabla", Tooltip -> Automatic,
    DisplayFunction :> (RowBox@{"(", #1, "\[CenterDot]\[Del]", ")", #2} &)]
MakeBoxes[Delta[a_], form:TraditionalForm] := TemplateBox[{MakeBoxes@Parenthesis@a}, "Delta", Tooltip -> Automatic,
    DisplayFunction :> (RowBox@{"\[Delta]", #1} &)]

And now for the most important part of the code, the function ExpandDerivative, which as the name suggest, expand the derivative-like functions:

(* Expand Derivatives/Vectors/Tensors on expr and apply custom rules *)
ExpandDerivative[expr_, rules_:{}] := expr //. Flatten@{
    (* Custom Rules *)
    rules,

    (* Linearity *)
    (op:grad|div|curl|Delta|Inactive[grad]|Inactive[div]|Inactive[curl]|Inactive[Delta]|Re|Im|Conjugate)[a_ + b__] :> op@a + op[+b],
    (op:Dp|Inactive[Dp]|Sum|Inactive[Sum])[a_ + b__, arg__] :> op[a, arg] + op[+b, arg],
    (op:Times|Dot|TensorProduct|Cross|DotNabla|Inactive[DotNabla])[a___, b_ + c__, d___] :> op[a, b, d] + op[a, +c, d],
    (op:grad|div|curl|Delta|Inactive[grad]|Inactive[div]|Inactive[curl]|Inactive[Delta]|Re|Im|Conjugate)[(op\[CapitalSigma]:Sum|Inactive[Sum])[a_, args__]] :> op\[CapitalSigma][op@a, args],

    (* Sum *)
    (op:Sum|Inactive[Sum])[s_*a_, v_Symbol] /; FreeQ[s, v] :> s*op[a, v],
    (op:Sum|Inactive[Sum])[s_, v_Symbol] /; FreeQ[s, v] :> s*op[1, v],

    (* Complexes *)
    Conjugate@(op:Times|Dot|Cross|TensorProduct)[a_, b__] :> op[Conjugate@a, Conjugate@op@b], (* Pass Conjugate to child *)
    (op:grad|div|curl|Delta|Inactive[grad]|Inactive[div]|Inactive[curl]|Inactive[Delta])[(opC:Re|Im|Conjugate)[a_]] :> opC@op@a, (* Pass Conjugate/Re/Im to parent *)
    Dp[(op:Re|Im|Conjugate)[a_], v_] :> op@Dp[a, v], (* Pass Conjugate/Re/Im to parent *)

    (* Triple Product *)
    Cross[a_?VectorQ, Cross[b_?VectorQ, c_?VectorQ]] :> b*a.c - c*a.b,
    Cross[Cross[a_?VectorQ, b_?VectorQ], c_?VectorQ] :> b*a.c - a*b.c,

    (* Quadruple Product *)
    Dot[Cross[a_?VectorQ, b_?VectorQ], Cross[c_?VectorQ, d_?VectorQ]] :> (a.c)*(b.d) - (a.d)*(b.c),

    (* Second Derivatives *)
    div@curl[_?VectorQ] :> 0,
    curl@grad[_?ScalarQ | _?VectorQ] :> 0,

    (* grad *)
    grad[(s_?ScalarQ) * (b_)] :> s*grad@b + b\[TensorProduct]grad@s,
    grad[(a_?VectorQ) . (b_?VectorQ)] :> a\[Cross]curl@b + b\[Cross]curl@a + DotNabla[a, b] + DotNabla[b, a], (* Use physics form *)
    grad[(s_?ScalarQ) ^ (n_?ConstantQ)] :> n*s^(n-1)*grad@s,
    grad[(n_?ConstantQ /; ScalarQ@n) ^ (s_?ScalarQ)] :> n^s*Log[n]*grad@s,

    (* div *)
    div[(s_?ScalarQ) * (b_?TensorQ)] :> s*div@b + b.grad@s,
    div[(a_?VectorQ) \[TensorProduct] (b_?VectorQ)] :> DotNabla[b, a] + a*div@b,
    div[(a_?VectorQ) \[Cross] (b_?VectorQ)] :> b.curl@a - a.curl@b,

    (* curl *)
    curl[(s_?ScalarQ) * (b_?VectorQ)] :> grad[s]\[Cross]b + s*curl@b,
    curl[(a_?VectorQ) \[Cross] (b_?VectorQ)] :> div[a\[TensorProduct]b - b\[TensorProduct]a],

    (* DotNabla *)
    DotNabla[(s_?ScalarQ) * (b_?VectorQ), c_?VectorQ] :> s*DotNabla[b, c],
    DotNabla[a_?VectorQ, (\[Beta]_?ScalarQ)*(c_?VectorQ)] :> c*a.grad@\[Beta] + \[Beta]*DotNabla[a, c],

    (* Dp *)
    Dp[(op:Times|Dot|Cross|TensorProduct)[a_, b__], v_Symbol] :> op[Dp[a, v], b] + op[a, Dp[op@b, v]],
    Dp[Power[a_?ScalarQ, b_?ScalarQ], v_Symbol] :> Power[a, b-1]*b*Dp[a, v] + Power[a,b]*Log[a]*Dp[b, v],

    (* Delta *)
    Delta[(op:Times|Dot|Cross|TensorProduct)[a_, b__]] :> op[Delta@a, b] + op[a, Delta@op@b],
    Delta@Power[a_?ScalarQ, b_?ScalarQ] :> Power[a, b-1]*b*Delta[a] + Power[a,b]*Log[a]*Delta[b]
}

Some examples. Calculating the divergent of Maxwell Stress Tensor in vaccum:

MST

Calculate the Einstein-Laub force density for linear dielectrics:

EL

Testing Poynting theorem in vaccum (no sources):

P

Where the first argument is the quantity being "tested".

Many other uses are possible and is fairly easy to extended some definitions.

POSTED BY: Thales Fernandes
Answer
4 months ago

enter image description here - Congratulations! This post is now a Staff Pick as distinguished on your profile! Thank you, keep it coming!

POSTED BY: Moderation Team
Answer
2 months ago

Very nice! But how does that relate to the built-in Symbolic Tensors functionality?

http://reference.wolfram.com/language/guide/SymbolicTensors.html

POSTED BY: Sam Carrettie
Answer
2 months ago

Mathematica built-in tensors are component based, whereas my simple implementation is coordinate free.

POSTED BY: Thales Fernandes
Answer
2 months ago

That is simply not true, both with regards to Mathematica and your own package.

Our symbolic tensors are very much coordinate and index free, and can be manipulated without ever introduce coordinates, indices, or what have you. What they lack is calculus support, which is a hard thing to do, even more so without a declarative syntax of the type you use. This is something we would like to circle back to, but there are many things we want to circle back to.

Given that you are defining Dp, which is an inherently coordinate-based (whether you define it as a non-tensorial local expression operator, or as a tensorial connection defined via a particular coordinate system), you clearly are not coordinate free. Where is your abstract connection, if you really want to be coordinate free?

POSTED BY: Itai Seggev
Answer
2 months ago

Our symbolic tensors are very much coordinate and index free, and can be manipulated without ever introduce coordinates, indices, or what have you.

My mind was on Calculus... Especifically Grad, Div, Laplacian, which are the tools I use and don't want specify clumsy coordinates to make the code unreadable...

Given that you are defining Dp, which is an inherently coordinate-based (whether you define it as a non-tensorial local expression operator, or as a tensorial connection defined via a particular coordinate system), you clearly are not coordinate free.

I use Dp not as coordinate derivative, but rather as parameter derivative. Specifically by assuming some deformation on the permittivity given by a parameter, which is the basis of Virtual Work.

Where is your abstract connection, if you really want to be coordinate free?

As I said, I called it "Package" with quotes for a reason. Why would I bother writing GR code if I don't do GR? My grad code is indeed coordinate free, for example. If I take the grad of the appropriate product of two tensors, it will give the correct output regardless of coordinate choice.

My main point was: give more power to Grad, Div, Laplacian, etc...

Mathematica Tensors may be coordinate free, but Mathematica Calculus definitely is not.

POSTED BY: Thales Fernandes
Answer
2 months ago

I see, so your interest was creating a vector calculus package for 3D Euclidean space, assuming you are always using an orthonormal basis. I agree this is a nice way to get there.

I would only add that there is a difference between "coordinate-free" and "index-free". Several of those packages you dismissed are in fact coordinate-free, even if they are not index-free. And the moment you have tensors of rank higer than 2, trying to stay index free is not so simple, or even particularly natural.

POSTED BY: Itai Seggev
Answer
2 months ago

I would only add that there is a difference between "coordinate-free" and "index-free".

Perhaps I used interchangeably component-free (or index-free) with coordinate-free... I should have made the distinction more clear.

I didn't use other packages because they are not component-free, and when using indices, things can get messier very quickly. And my code is just Leibniz rule for tensors, there is no need for any component to appear.

And I didn't use Mathematica own vector derivatives because they are not coordinate-free (in the sense I need to specify the coordinate, obviously I can derive in any coordinate with them, which is pretty handy). It is possible to extend this "Package" for more definitions, for example, Riemann curvature tensor:

R

Which is component-free and coordinate free. But it would not be too useful since the maximum you can do without specifying components or defining derivative-like operators is Leibniz rule.

Note: The only book I have found so far with good component-free notation in GR is: Grøn, Hervik - Einstein's General Theory of Relativity With Modern Applications in Cosmology. The math is more complicated and sometimes is necessary to revert back to index notation.

POSTED BY: Thales Fernandes
Answer
2 months ago

You might want to take a look at the Atlas2 add-on:

http://www.digi-area.com/Mathematica/atlas/

POSTED BY: Murray Eisenberg
Answer
2 months ago

Group Abstract Group Abstract