Exterior Differential Calculus
and
Symbolic Matrix Algebra

version 3.6.0 - May 2007

S. Bonanos

Institute of Nuclear Physics,  NCSR "Demokritos", 
Aghia Paraskevi 15310, Greece. 

Please send comments, suggestions, bug reports to:  sbonano@inp.demokritos.gr

Overview

This package enables Mathematica to carry out calculations with differential forms. It defines the two basic operations - Exterior Product (Wedge) and Exterior Derivative (d) - in such a way that:
 (1)  they can act on any valid Mathematica expression
 (2)  they allow the use of any symbols to denote differential forms
 (3)  input - output notation is as close as possible to standard usage

There are two versions of this package: scalarEDCcode.nb and matrixEDCcode.nb. The first can handle scalar differential-form expressions only, while the second can also handle matrix-valued differential forms, i.e., matrices whose components are (scalar) differential forms.  matrixEDCcode.nb can be used for calculations with scalar forms, but it is  ~15% slower than  scalarEDCcode.nb.
 	The matrix package, allowing user-controlled application of trace identities and the Cayley-Hamilton theorem, can also be used for symbolic matrix calculations (see sample calculations (2) and (3) in the last section).

For exterior multiplication we use the infix operator symbol ⋀ (Wedge). This symbol is included in Mathematica (see Sections 3.10.1, 3.10.4) but does not have a pre-defined meaning. It can be input by typing  Esc ^ Esc  (where ^ is the symbol for exponent) or  \ [Wedge] (without space after \). Finally, it can be input by clicking the appropriate button in the notebook SymbolPalette.nb. 
  	The same symbol (instead of Dot) is used to denote matrix multiplication.

All expressions are treated as sums of terms, which may have factors that are special symbols (declared by the user as denoting differential forms or matrices) or have Head Wedge or d, and rules are given for manipulating such expressions. All other (undeclared) symbols are treated as scalar 0-forms. The ordering of terms in an expression is determined by the built-in Mathematica rules for Plus and Times, which sometimes results in an unexpected ordering of terms. All terms in a sum, and all components of a matrix, must be of the same (differential-form) degree.

Warning: The executable cells in this Manual must be evaluated in the order they appear, after evaluating matrixEDCcode.nb first.

Important Note: The notebooks are in the form of "open code", and the many global variables introduced are not protected. This is intentional, in order to allow the user to add new definitions on the fly, as required in a calculation, without having to go through Unprotect. However, it also carries the danger that an unprotected symbol might be used in a way that was not intended. If one obtains results that don't make sense, it is most likely that some global variable is used inappropriately (e.g.,  a symbol, declared to represent a matrix, is used in a later calculation as a scalar).  Reevaluating the relevant EDCcode.nb notebook clears all variables and allows a fresh start.

Acknowledgement: I wish to thank Wolfram Research for a Visiting Scholar Grant, during which version 2.6.1 of this package was developed.

The present version - 3.6.0 - is an update to make EDC compatible with Mathematica 6. See version history in the EDC web page.

scalarEDC

The notebook scalarEDCcode.nb allows calculations with scalar differential forms only. The following functions / operations are defined (examples are given in the cells in brackets. Sections marked by  ■  can be skipped on a first reading):

Wedge[x_,y_,...]  or   x_⋀y_⋀... :: gives the exterior product of the expressions x,y,... Wedge reduces to Times when acting on 0-forms.

Wedge[2  a d[x] + 3  b d[y] + 5  c d[z], Cos[a x + b y] d[y]]

Note that Wedge products are ordered alphabetically.  Instead of Wedge[...] one can use the infix symbol ⋀:

(2  a d[x] + 3  b d[y] + 5  c d[z])⋀(Cos[a x + b y] d[y])

d[x_] :: gives the exterior derivative of the expression x, assuming all symbols are variable.

d/@{A, A[0], A[x], A[2][x, y]}

d[%%]

Constants must be explicitly declared:

d[a] = 0; d[b] = 0; %

d[x^n]

d[n] = 0; %

reWrite[x_] :: rewrites the expression x, collecting together all terms having identical Wedge[_] or d[_] factors, applies the simplification rules defined in the list simpRules, and factors the coefficient of each such Wedge[_] or d[_] term. In complicated expressions, one can define a set of variables in the list varList to be used to collect terms in these coefficients so as to further simplify the result - see next section. The proper choice of rules in simpRules and of variables in varList is crucial for the effective simplification of complicated expressions. The default values of the global variables simpRules and varList are: simpRules = {x_. Cos[z_]^2 + x_. Sin[z_]^2 -> x}; varList = {};

reWrite[%%% - 3  %%%*a/5]

Another application of reWrite, where the default value of simpRules is applied:

d[x]⋀d[y]⋀d[z] /. {x -> r Sin[θ] Cos[ϕ], 
        y -> r Sin[θ] Sin[ϕ], z -> r Cos[θ]}

reWrite[%]

The use of varList

Suppose the function  g  of x, y, z  is known to depend on x, y through the two particular combinations - 
(x + y)/((x^2) + (y^2)) and (x y + 1)/((x^2) + (y^2)) -  and g is written with these as arguments. Then

dg = 
      d[g[(x + y)/((x^2) + (y^2)), (x y  + 1)/((x^2) + (y^2)), 
          z]]

reWrite[dg]

Note that reWrite collects together all terms having  d[x], d[y] as factors, but it does not group together terms with the same derivatives of g in the coefficients of d[x], d[y]. Therefore, choose

varList = {Derivative[__][g][__]};

Now reWrite gives:

reWrite[dg]

Verify that it is a closed form:

reWrite[d[%]]

varList = {};

FormDegree[x_]:: returns the (differential-form) degree of the expression x. It assumes that all terms in a sum have the same degree.

DeclareForms[{p1_,x__},{p2_,y__},...] :: allows the user to introduce his particular notation for fundamental (undefined) p-forms. p1, p2 must be positive integers (denoting the degree of the differential form(s) x__, y__) and x__, y__  sequences of one or more symbols separated by commas. For example, evaluating the command DeclareForms[{1,a1,q[_]},{2,a2,b2,c2}] implies that, in all subsequent calculations, the symbols a1, q[1], q[2],  etc., will be treated as 1-forms and the symbols a2, b2, c2  as 2-forms. If only forms of the same degree are introduced, the curly brackets can be omitted: DeclareForms[1,a1,q[_]] is the same as DeclareForms[{1,a1,q[_]}]. DeclareForms also assigns the appropriate FormDegree to the newly introduced forms, and makes appropriate changes to other internal variables. All other undefined symbols are treated as scalars (0-forms). 
Note: Only undefined symbols should be declared as forms using DeclareForms. If all forms are expressed in terms of the coordinate differentials (d[x], d[x]⋀d[y], etc.) DeclareForms should not be used. The effects of DeclareForms  cannot be undone.

DeclareForms[{1, a1, q[_]}, {2, a2, b2, c2}]

q[1]⋀a1

Wedge orders forms in canonical order by degree (lower order forms first):

c2⋀%⋀a2

One can give arbitrary, but consistent, definitions for the exterior derivatives of undefined forms:

d[q[1]] = q[2]⋀q[3]; d[q[2]] = q[3]⋀q[1]; 
    d[q[3]] = q[1]⋀q[2];

d[{q[2]⋀q[3], q[3]⋀q[1], q[1]⋀q[2]}]

Warning: One must not use Set (=) for relations between udefined forms. For example, if one discovers, in the course of a calculation, that the equations satisfied by an undefined form (say q[3]) imply that it can be expressed in terms of scalars as, say,  u d[v] + d[w], one must avoid making the assignment q[3] = u d[v] + d[w] because it may lead to infinite loops (depending on previous calculations using q[3]). One must instead use the replacement rule q[3] -> u d[v] + d[w] where needed, or restart the calculations using a different notation for q[3] (say q3, which must not be declared as an 1-form), and then setting  q3 = u d[v]+d[w].

Bar[x_]:: returns the complex conjugate of the expression x, assuming all symbols are complex. If some symbols are real, the user must set Bar[ai]=ai, for each real symbol ai. If most symbols are real one can define Bar[x_Symbol]:=x, and still set Bar[ci]=ciB;Bar[ciB]=ci; for pairs of complex-conjugate symbols ci,ciB.

r1 = {z -> x + I y}

r2 = Join[r1, Bar[r1]]

Bar[x_Symbol] := x; Bar[z] = zB; Bar[zB] = z; r2

ScalFormQ[x_], BasicScalFormQ[x_]: return True if x is a (basic) differential-form expression and False otherwise. Basic means not expressible (as sum or product) in terms of other forms.

testExpr = {x, d[x], 3 d[x], a1, a1 + q[1], a1⋀q[1]};

\(\(ScalFormQ/@testExpr
\)\)

\(\(BasicScalFormQ/@testExpr
\)\)

Calculations with Series

When working with series expansions, one must keep in mind that the notation d[x] for the differential of the expansion symbol (coordinate) cannot be used because Series will consider d[x] a function of x and formally expand it (using d[0]=0):

Series[d[x]/(1 - x), {x, 0, 3}]

One must therefore use DeclareForms to introduce a new symbol for the differential of that scalar variable (and appropriately define its exterior derivative). For example:

\(DeclareForms[1, dx]; 
d[x] = dx; d[dx] = 0; \)

Then one can do calculations with series (one can still input d[x]):

ser1 = Series[d[x]/(1 - x), {x, 0, 3}]

d[ser1]

ser1⋀ser1

d[y^2 ser1]

Series has been modified to automatically expand (to the appropriate order) series coefficients that depend on x:

% /. y -> Sqrt[1 - x] u

matrixEDC

The notebook matrixEDCcode.nb allows calculations with differential form matrices - both explicit (list of lists) and symbolic (special symbols denoting matrices). In addition to the functions defined in scalarEDCcode.nb, which can now have matrices as arguments (except for DeclareForms), the following functions / operations are defined (examples are given in the cells in brackets. Sections marked by  ■  can be skipped on a first reading):

DeclareMatrixForms[{p1_,x_,xt_},{p2_,y_,yt_},...] :: allows the user to introduce his particular notation for fundamental (undefined) p-form matrices. p1 must be a positive integer or zero (denoting the differential form degree of x_, xt_) and x_, xt_  are symbols denoting a p1-form matrix and its transpose, respectively. For example, evaluating the command DeclareMatrixForms[{0,S,S},{1,U,-U},{1,Q,Qt}] has the effect that, in all subsequent calculations, the symbol  S  will be treated as a 0-form symmetric matrix, the symbol  U  as an 1-form antisymmetric matrix, and the symbols  Q, Qt  as 1-form matrices which are transposes of each other.

DeclareMatrixForms[{0, S, S},{1,U,-U},{1, Q, Qt}]

Transpose[S⋀U]

FormDegree[%]

One can also use indexed variables to denote sets of matrices with common properties. For example, after evaluating the command

DeclareMatrixForms[{0, S0[_], S0[_]}, {1, A1[_], -A1[_]}, {1, 
          N1[_], N1t[_]}]

the symbols S0[1], S0[2],... are treated as symmetric 0-form matrices, the symbols A1[1], A1[2],... as antisymmetric 1-form matrices, and the symbols N1[1], N1[2],..., N1t[1], N1t[2],..., as 1-form matrices which are, respectively, transposes of each other.

Transpose[N1[2]⋀S0[1]⋀N1[1]]

The minus sign is correct - the N1[i] are 1-form matrices!
Note : No special symbol is used for the identity matrix, as it is not needed in symbolic calculations. The functions tr and toComponents have been defined so as to take into account the missing identity matrix.

Transpose[x_] :: the built-in definition of Transpose is extended so that it handles properly symbolic matrix expressions containing the symbols introduced with DeclareMatrixForms.

SymQ[x_] :: returns True if Transpose[x]=x and False otherwise.
AntisymQ[x_] ::returns True if Transpose[x]=-x and False otherwise.

In SymQ/AntisymQ  equality is determined after taking into account the simplification rules contained in simpRules.

Transpose[
      (Cos[a]^2 + Sin[a]^2) S0[1]⋀A1[1] + A1[1]⋀S0[1]]

AntisymQ[%]

MatQ[x_], BasicMatQ[x_]: return True if x is a (basic) matrix expression and False otherwise. Basic means not expressible (as sum or product) in terms of other matrices.

testExpr
       = {x, d[x], S, d[S], N1[2] + N1t[2], S⋀U, d[x]⋀U, 
        Bar[N1[_]]};

\(\(MatQ/@testExpr
\)\)

\(\(BasicMatQ/@testExpr
\)\)

tr[x_] :: returns the trace of the matrix expression x. If  x is a scalar, it is assumed to multiply the identity matrix, and thus tr returns x*matDimension, where matDimension is a global variable to be assigned by the user. 
When tr  acts on products of symbolic matrices, it can apply the cyclic property of trace (as well as the fact that tr[Transpose[x]]=tr[x]), to reduce its argument to a "canonical" form. It can thus be made to recognize that the trace of the matrix expression 2 A⋀B⋀C-C⋀A⋀B-At⋀Ct⋀Bt is identically zero (for 0-form matrices A,B,C with transposes At,Bt,Ct respectively).

DeclareMatrixForms[{0, A, At}, {0, B, Bt}, {0, C, Ct}]

The behavior of tr is controlled by the value of the global variable traceOption,  as can be seen in the following examples:

testExpr={At,Ct⋀Bt,At⋀Ct⋀Bt,At⋀Ct⋀Bt⋀At,2 A⋀B⋀C-C⋀A⋀B-At⋀Ct⋀Bt};

traceOption = 0; tr/@testExpr(** No rules applied, other than linearity **)

traceOption = 1; tr/@testExpr(** Only cyclic rule applied **)

traceOption = 2; tr/@testExpr(** Cyclic and Transpose rules applied **)

Note that, for single matrices, tr[At]=tr[A] independently of traceOption. 
	Warning: The default value of traceOption  is 2, to allow maximum simplification of complicated  tr  expressions. However, this value can sometimes lead to infinite loops if used together with a commutator rule like S0[1]⋀S0[3]->S0[3]⋀S0[1]+2I S0[2] which reverses the canonical order of matrix products.

tr applies other rules also; for example, since S is a symmetric 0-form matrix and U antisymmetric,

tr/@{S⋀U,S⋀U⋀S⋀S}

toComponents[x_,rul_List]:: returns the component form of the symbolic matrix expression x, introducing the identity matrix where needed. The components of the basic symbolic matrices must be given as a set of replacement rules in the list rul and must satisfy the assumed matrix symmetries. If x is a scalar, toComponents returns x*IdentityMatrix[matDimension].
Note : If the components of a p-form matrix are undefined scalar p-forms, DeclareForms must be used to introduce them. For example, if U has been declared an antisymmetric 1-form matrix, and rul gives the components of U  in terms of the undefined scalar 1-forms u1, u2, u3, one must evaluate DeclareForms[1,u1,u2,u3] before using toComponents.

matDimension = 3;
matComp={A -> Array[Ac, {3, 3}],At -> Transpose[Array[Ac, {3, 3}]],S -> Array[Sc, {3, 3}],U->Array[ucomp, {3, 3}]}

The components must be made to satisfy the assumed matrix symmetries, and must be declared to have the appropriate differential form degree (if >0):

\(Sc[i_, j_] := Sc[j, i] /; i > j; 

    ucomp[i_, j_] := Sum[Signature[{i, j, k}] uc[k], {k, 3}]\)

DeclareForms[1, uc[_]]

matComp

MatrixForm[toComponents[S⋀U,matComp]]

tr[%]

FormDegree[%%]

toComponents understands that, in the expression  A + lamda, the scalar lamda multiplies the identity matrix:

MatrixForm[toComponents[A  +  lamda, matComp]]

tr[%]

(A  +  lamda)⋀Transpose[A  +  lamda]

tr[%]

MatrixForm[toComponents[%%, matComp]]

SymQ[%]

tr[%%]

Collect[%, lamda, Factor]

matDimension =.

CHrul[n_Integer, mat_:$]:: returns a rule for implementing the Cayley-Hamilton theorem for the n x n  0-form matrix mat. The second argument is optional; if it is missing CHrul returns a delayed rule that can be applied to any n x n  0-form matrix. CHrul introduces the symbol $ as a 0-form matrix. After calling CHrul, matDimension is set to n.

CHrul[3, A]

CH3 = CHrul[3]

The rule applies only when $ is a 0-form matrix expression. The pattern matches products of matrices also:

reWrite[A⋀B⋀A⋀B⋀A⋀B⋀A /. CH3]

CH4 = CHrul[4]

Now matDimension is set to 4, so CH3 has no effect:

reWrite[A⋀B⋀A⋀B⋀A⋀B⋀A /. CH3]

Warning : The rules returned by CHrul must never be applied repeatedly (//.), as there is a term on the rhs that matches the lhs (tr[lhs]),  To apply these rules repeatedly, one must use another symbol, say trnA, for the trace of the product of n A matrices.

TimeConstrained[A⋀A⋀A //. CHrul[3, A], 3]

Timing[
      A⋀A⋀A //. 
        (CHrul[3, A] /. tr[A⋀A⋀A] -> tr3A)]

Note : Wedge, instead of Dot,  must be used to multiply matrices together and scalar differential forms with matrices. Wedge can handle both symbolic and explicit matrices (list of lists). Wedge also gives the correct result when multiplying explicit matrices (two-dim arrays) with vectors (one-dim arrays), but operations with symbolic vectors are not supported. For example:

Mat = A /. matComp

Vec = Array[v, {3}]

DeclareForms[1, v[_]]

(i.e., Vec is an 1-form vector)

Mat⋀Vec

Vec⋀Mat

reWrite[Vec⋀Mat⋀Vec]

Vec⋀d[Vec]

Vec⋀Mat⋀d[Vec]

d[x]⋀Vec

d[Mat]⋀Vec

reWrite[Vec⋀%]

Introducing Custom Notation : Sometimes one may want to use a shorthand notation of the form K[A] for an operation with matrices (or other variables) that should not be interpreted as a function under d. I.e., one wants to avoid the following

d[{Inverse[S], MyNotation[form, vars]}]

This can be achieved by: 	(1) Calling the function NoDif[heads__], 
				(2) Setting the value of FormDegree[..] if > 0 and
				(3) Setting the values of BasicMatQ[..] / BasicScalFormQ[..] if True. 
For example,

\(NoDif[Inverse, MyNotation]; 
BasicMatQ[Inverse[x_]] := BasicMatQ[x]; 

    FormDegree[MyNotation[x_, y_]] := FormDegree[y]; \)

Now d has no effect

d[{Inverse[S], MyNotation[form, vars]}]

One can then assign values for the exterior derivatives of these operations and make other definitions:

\(d[Inverse[x_]] := 
      \(-\((Inverse[x]⋀d[x]⋀Inverse[x])\)\) /; 
        BasicMatQ[x] && FormDegree[x] === 0; 

    x_⋀Inverse[x_] := 1 /; BasicMatQ[x] && FormDegree[x] === 0; 
    Inverse[x_]⋀x_ := 1 /; BasicMatQ[x] && FormDegree[x] === 0; 

    S⋀d[Inverse[S]]\)

Of course, for a single matrix, it is much simpler to use a new symbol  DeclareMatrixForms[{0,invS,invS}] and then set d[invS]=-(invS⋀d[S]⋀invS);invS⋀S=1;S⋀invS=1;!

Sample Applications in Mathematical Physics

Important: Before proceeding, evaluate matEDCcode.nb to clear all previously defined symbols.
	(All timings refer to an Apple iBook G3/900 computer running Mathematica 3.0 under Mac OS 9.2.2)

(1) Computing all Newman-Penrose scalars for the Vaidya metric (~ 25 sec)

The following code calculates all Newman-Penrose scalars for the Vaidya metric and checks the cyclic and Bianchi idendities.
(Takes  ~ 25 sec with matEDCcode and ~ 21 sec with scalEDCcode. More than half of this time is spent in checking the identities). 
To do the same calculation for other metrics make appropriate changes to the null-frame 𝕝,𝕞,𝕞B,𝕟 definitions. For complicated metrics additional simplification statements (or rules in simpRules) may be needed.

Off[General::"spell", General::"spell1"]; Bar[x_] := x; 
    sT = TimeUsed[];

DEFINE NULL FRAME  𝕝,𝕞,𝕞B,𝕟

Delta0 = (r^2) + (a^2) ; Delta = Delta0 - 2  M r; 
      Sigma = (r^2) + (a^2) (x^2); (** x=Cos[θ] **)

BaseFrame = 
      reWrite[{d[t] - a (1 - (x^2)) d[φ], Sigma/Delta d[r], 
          -Sigma d[x]/Sqrt[1 - (x^2)], 
          a d[t] - Delta0 d[φ]}];

(** Transform t, φ coordinates **)
d[t] = d[u] + Delta0 d[r]/Delta; 
      d[φ] = d[nφ] + a d[r]/Delta;

\(lCoef = {1, \(-1\), 0, 0}; 

    𝕝 = reWrite[lCoef.BaseFrame]\)

\(nCoef = {1, 1, 0, 0} Delta/\((2\ Sigma)\); 

    𝕟 = reWrite[nCoef.BaseFrame]\)

\(mCoef = {0, 0, \(-1\), I\ Sqrt[1 - x^2]}/Sqrt[2\ Sigma]; 

    𝕞 = reWrite[mCoef.BaseFrame]\)

𝕞B = Bar[𝕞];

LineElement = 
      reWrite[2  𝕝*𝕟 - 
          2  𝕞*𝕞B]

d[a] = 0; M = m[u]; (** Vaidya - for Kerr set d[M]=0 **)

END DEFINE METRIC

iden = 
      reWrite[I 
          𝕝⋀𝕟⋀𝕞⋀𝕞B]

star4 = Times[u__, Wedge[x_]] :> Times[u]

vol = iden /. star4

Define the self-dual 2-forms f[i]:

f[0] = reWrite[𝕞⋀𝕝]

f[1] = 
      reWrite[(
            𝕞⋀𝕞B - 
              𝕝⋀𝕟)/2]

f[2] = reWrite[𝕟⋀𝕞B]

fB[i_] := fB[i] = Bar[f[i]];

df[i_] := df[i] = reWrite[d[f[i]]];

CALCULATE THE SPIN COEFFICIENTS:

NPκ = 
      -reWrite[I 𝕝⋀df[0]/vol] /. star4

NPσ = 
      -reWrite[I 𝕞⋀df[0]/vol] /. star4

NPρ = 
      -reWrite[I 𝕝⋀df[1]/vol] /. star4

NPμ = 
      -reWrite[I 𝕟⋀df[1]/vol] /. star4

NPπ = 
      -reWrite[I 𝕞B⋀df[1]/vol] /. star4

NPτ = 
      -reWrite[I 𝕞⋀df[1]/vol] /. star4

NPν = 
      -reWrite[I 𝕟⋀df[2]/vol] /. star4

NPλ = 
      -reWrite[I 𝕞B⋀df[2]/vol] /. star4

NPγ = 
      Factor[(
            NPμ - reWrite[I 𝕞⋀df[2]/vol] /. star4)
            /2]

NPϵ = 
      Factor[(
            NPρ - reWrite[I 𝕞B⋀df[0]/vol] /. 
              star4)/2]

NPβ = 
      Factor[(
            NPτ - reWrite[I 𝕟⋀df[0]/vol] /. star4)
            /2]

NPα = 
      Factor[(
            NPπ - reWrite[I 𝕝⋀df[2]/vol] /. star4)
            /2]

Definition of connection 1-forms:

ω0 = 
      reWrite[NPτ*𝕝 + NPκ*𝕟  - 
          NPρ*𝕞 - NPσ*𝕞B]

ω1 = 
      reWrite[NPγ*𝕝 + NPϵ*𝕟  - 
          NPα*𝕞 - NPβ*𝕞B]

ω2 = 
      reWrite[NPν*𝕝 + NPπ*𝕟  - 
          NPλ*𝕞 - NPμ*𝕞B]

TimeUsed[] - sT

Definition of curvature 2-forms:

Ω[0] = d[ω0] + 2 ω0⋀ω1;

Ω[1] = d[ω1] + ω0⋀ω2;

Ω[2] = d[ω2] + 2 ω1⋀ω2;

Do[Ω[i] = reWrite[Ω[i]], {i, 0, 2}]

TimeUsed[] - sT

WEYL COMPONENTS AND SYMMETRIES

Ψ0 = 
      -reWrite[I f[0]⋀Ω[0]/vol] /. star4

Ψ1 = 
      -reWrite[I f[1]⋀Ω[0]/vol] /. star4

reWrite[
      Ψ1 + reWrite[I f[0]⋀Ω[1]/vol] /. 
        star4]

Ψ3 = 
      -reWrite[I f[2]⋀Ω[1]/vol] /. star4

reWrite[
      Ψ3 + reWrite[I f[1]⋀Ω[2]/vol] /. 
        star4]

Ψ4 = 
      -reWrite[I f[2]⋀Ω[2]/vol] /. star4

fΩ02 = 
      reWrite[I f[0]⋀Ω[2]/vol] /. star4;

fΩ20 = 
      reWrite[I f[2]⋀Ω[0]/vol] /. star4;

fΩ11 = 
      reWrite[I f[1]⋀Ω[1]/vol] /. star4;

Ψ2 = 
      -reWrite[
          (4  fΩ11 + fΩ02 + fΩ20)
              /6]

reWrite[fΩ02 - fΩ20]

NON-ZERO RICCI COMPONENTS

Λ = 
      reWrite[(
            2  fΩ11 - fΩ02 - fΩ20)
            /6]

Φ[i_, j_] := 
      Φ[i, j] = 
        reWrite[I Ω[i]⋀fB[j]/vol] /. star4;

cnt = 0; 
    Do[tmp = Φ[i, j]; 
      If[tmp =!= 0, 
        Print["Φ[", i, ",", j, "] = ", tmp]; 
        cnt++], {j, 0, 2}, {i, 0, 2}]; 
    If[cnt === 0, 
      If[Λ === 0, Print["Ricci Flat"], 
        Print["Einstein Space"]]];

TimeUsed[] - sT

VERIFY Φ SYMMETRY

reWrite[Φ[0, 1] - Bar[Φ[1, 0]]]

reWrite[Φ[0, 2] - Bar[Φ[2, 0]]]

reWrite[Φ[1, 2] - Bar[Φ[2, 1]]]

reWrite[Φ[0, 0] - Bar[Φ[0, 0]]]

reWrite[Φ[1, 1] - Bar[Φ[1, 1]]]

reWrite[Φ[2, 2] - Bar[Φ[2, 2]]]

TimeUsed[] - sT

VERIFY CYCLIC IDENTITY

reWrite[
      (Ω[1] - Bar[Ω[1]])
            ⋀𝕞 + 
        Bar[Ω[2]]⋀𝕝 - 
        Ω[0]⋀𝕟]

reWrite[
      -(Ω[1] + Bar[Ω[1]])
              ⋀𝕝 + 
        Bar[Ω[0]]⋀𝕞 + 
        Ω[0]⋀𝕞B]

reWrite[
      -(Ω[1] + Bar[Ω[1]])
              ⋀𝕟 + 
        Bar[Ω[2]]⋀𝕞B + 
        Ω[2]⋀𝕞]

TimeUsed[] - sT

VERIFY BIANCHI IDENTITIES

reWrite[
      d[Ω[0]] + 
        2 (ω0⋀Ω[1] - 
              Ω[0]⋀ω1)]

reWrite[
      d[Ω[1]] + ω0⋀Ω[2] - 
        Ω[0]⋀ω2]

reWrite[
        d[Ω[2]] - 
          2 (ω2⋀Ω[1] - 
                Ω[2]⋀ω1)] // Timing

TimeUsed[] - sT

On[General::"spell", General::"spell1"];

More examples of calculations in Riemannian Geometry can be found in the package Riemannian Geometry and Tensor Calculus (Sections 5 and 6).

(2) Polynomial Matrix Algebra

Suppose one wishes to simplify products of polynomials in a matrix variable X, say in 3 dimensions, using the Cayley-Hamilton theorem. matrixEDCcode cannot easily handle the notation X^n for powers of a matrix, but one can use the notation X[n] for the nth power of the matrix X[1] and define appropriate rules for wedge-products of  X[k_] matrices. Thus

DeclareMatrixForms[{0, X[_], Xt[_]}]

X[i_]⋀X[k_] := X[i + k];

Use the Cayley-Hamilton Theorem to reduce the X[k]

CH3 = CHrul[3, X[1]]

CH3rhs = CH3[[1]][[2]]

For simplicity (and to avoid infinite loops!) set

tr[X[1]] = s1; tr[X[2]] = s2; tr[X[3]] = s3;

CH3rhs

In this notation X[0] is the identity matrix. One can either set X[0]=1, or  set

tr[X[0]] = 3; Unprotect[Transpose]; Transpose[X[0]] = X[0];  
    Protect[Transpose];

and use it to rewrite CH3rhs:

CH3rhs = reWrite[CH3rhs⋀X[0]]

reWrite[tr[%]](** check that trace of CHrhs = tr[X[3]] **)

Now define rules for simplifying X[i] matrices

X[i_ /; i > 2] := reWrite[X[i - 3]⋀CH3rhs]

With these rules one can calculate arbitrary powers of matrix polynomials:

X[10] // Timing

reWrite[X[5]⋀X[7]] // Timing

matPoly = X[0] + k1 X[1] + k2 X[2]

reWrite[matPoly⋀matPoly⋀matPoly] // Timing

Finding the inverse of  matPoly

testPoly = Sum[a[i] X[i], {i, 0, 2}]

zer0 = reWrite[testPoly⋀matPoly - X[0]]

eqs = {zer0 /. X[0] -> 1, zer0 /. X[1] -> 1, zer0 /. X[2] -> 1} /. 
        X[_] -> 0

invPolyRul = Solve[eqs == 0, {a[0], a[1], a[2]}][[1]]

invPolyRul = reWrite[invPolyRul]

invPoly = reWrite[testPoly /. invPolyRul]

reWrite[matPoly⋀invPoly](** check inverse **)

(3)  SU(2) algebra: Trace calculations

Off[General::"spell", General::"spell1"];

Introduce notation and define properties of Pauli matrices:

DeclareMatrixForms[{0, σ[_], σt[_]}]

tr[σ[a_]] = 0;

matDimension = 2;

\(delta[a_, b_] := If[a === b, 1, 0]; 

    epsilon[a_, b_, c_] := Signature[{a, b, c}]; \)

σ[a_]⋀σ[b_] := 
      delta[a, b] + I Sum[epsilon[a, b, c] σ[c], {c, 3}]

Define A, B  to be two SU(2) algebra vectors:

A = Sum[a[i] σ[i], {i, 3}]

B = Sum[b[i] σ[i], {i, 3}]

reWrite[A⋀B]

reWrite[tr[A⋀B]]

reWrite[tr[A⋀A⋀B]]

reWrite[tr[A⋀A⋀B⋀B]]

reWrite[tr[A⋀B⋀A⋀B]]

reWrite[tr[A⋀A⋀B⋀B⋀A⋀B]] // Timing

matDimension =. ; On[General::"spell", General::"spell1"];