Exterior Differential Calculus
and
Symbolic Matrix Algebra

version 3.8.9 - May 2013

S. Bonanos

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

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

1  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

Another use of this package is for doing algebraic / differential calculations with "symbolic matrices", i.e., with symbols satisfying special multiplication rules, which can be interpreted as representing matrices, quantum operators, Lie algebra generators, Maurer-Cartan forms etc. Any symbol can be defined to be a "symbolic matrix", i.e., to have special multiplication properties. But in this case the user must give the extra multiplication (Wedge) rules that define his/her problem (see Examples).

The infix operator symbol ⋀ (Wedge) has been defined as a generalized associative, non-commutative multiplication operator. Its action depends on the properties of the symbols it acts on: it becomes exterior multiplication on differential forms, matrix multipication on (symbolic or explicit) matrices and reduces to Times for other symbols.  The symbol ⋀ (Wedge) is included in Mathematica  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 palette EDCpalette.nb.

The letter "d" has been defined to act as the exterior differential operator. Thus one can use the standard notation  d[x], d[y], d[z], d[x]⋀d[y], d[x]⋀d[y]⋀d[z], d[f[x,y]] for differential forms and for exterior differentiation.

All expressions (or list entries) 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.

Important Note: The new functions (Wedge, d, FormDegree, Bar, MatQ, etc.) and global variables (simpRules, varList, etc.) introduced in the package 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 initialization code matrixEDC.m clears all variables and allows a fresh start.

In previous releases there were two versions of this package: scalarEDCcode.nb and matrixEDCcode.nb. As the efficiency of EDC code has improved (and personal computers have become much faster!), the only advantage of the scalar package (15% faster execution) no longer justifies maintaining this version. Thus, the new version is only the matrix version in package form (matrixEDC.m), which, of course, can be used for scalar calculations also.

The executable cells in this Manual must be evaluated in the order they appear, after placing the file matrixEDC.m in one of the directories that  Mathematica looks at when searching for extra packages. Such directories can be found by evaluating the command $Path:

$Path

Load initialization code:

<< matrixEDC.m

$EDCversion

2  EDC Function Definitions

2.1  The three basic operations: Wedge, d, reWrite

Let us define two differential forms in terms of the coordinate differentials in 3-dimensional Euclidean space:

f1 = 2  a d[x] + 3  b d[y] + 5  c d[z];

f2 = a Cos[a x + b y] d[x] - b Sin[a x + b y] d[y];

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

f12 = Wedge[f1, f2]

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

f1⋀f2

•	The function reWrite[x_] collects identical Wedge[__] or d[_] terms together and factors the coefficients of such terms:

reWrite[f12]

•	The function 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], u[x, y]/v[x], Cos[a x]}

Constants must be explicitly set to have 0 exterior derivative:

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

d[f2]

df2 = reWrite[%]

reWrite[d[f1]⋀f2 - f1⋀d[f2] - d[f12]]

•	Wedge, d, and reWrite can also be used with matrices; Wedge on matrices gives matrix multiplication:

Let us define the rotation matrix:

R3 = {{Cos[a x + b y], -Sin[a x + b y], 0}, {Sin[a x + b y], 
          Cos[a x + b y], 0}, {0, 0, 1}}

Wedge[R3, R3]

d[R3]

Wedge can be used to multiply a scalar form with a matrix:

reWrite[d[x]⋀%]

Wedge[Inverse[R3], d[R3]]

reWrite[%]

•	Wedge can take any number of arguments; it reduces to Times when fewer than 2 of its arguments are differential forms or matrices:

Wedge[p, q^2, r^3, f1]

Wedge[p, q^2, r^3, R3]

2.2  Defining symbolic forms and matrices: DeclareForms and DeclareMatrixForms

In addition to explicit forms and matrices (like f1, f2, R3 above), EDC allows the user to introduce symbolic forms and matrices (of arbitrary differential-form degree); Wedge and d will then act correctly on these symbols. Any symbol can be defined to be a symbolic form or matrix using the commands DeclareForms and DeclareMatrixForms.

(i) DeclareForms

The command DeclareForms[{p1_,x__},{p2_,y__},...]  allows the user to introduce his particular notation for fundamental (=not expressible in terms of other forms) p-forms. DeclareForms takes any number of arguments of the form  {p1_,x__}, where p1 is a positive integer (denoting the degree of the differential form(s) x__) and x__  a sequence 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[_]}].

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 (alphabetical) order by degree (lower order forms first):

c2⋀%⋀a2

d  takes into account the degree of the symbolic forms:

d[a1⋀q[1]]

d[a2⋀b2]

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].

(ii) DeclareMatrixForms

The command DeclareMatrixForms[{p1_,x_,xt_},{p2_,y_,yt_},...] allows the user to introduce his particular notation for fundamental (undefined in terms of other matrices) symbolic p-form matrices. DeclareMatrixForms takes any number of arguments of the form  {p1_,x_,xt_}, where p1 is 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,W,Wt}] 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  W, Wt  as 1-form matrices which are transposes of each other.

Note: Only undefined symbols should be declared as matrices using DeclareMatrixForms. If one introduces a new symbol for a matrix that is defined in terms of other matrices (say, P=S⋀U, or explicitly in terms of a list of components), that symbol (P) must not be declared a symbolic matrix using DeclareMatrixForms. The effects of DeclareMatrixForms  cannot be undone.

DeclareMatrixForms[{0, S, S},{1,U,-U},{1, W, Wt}]

The built-in function Transpose  (see section 2.4) has been extended to apply elementary identities to symbolic matrices, e.g.,  Transpose[A⋀B]=(-1)^(p*q)Transpose[B]⋀Transpose[A], where A, B are p, q-form matrices, respectively:

Transpose[S⋀U]

d[%]

Transpose[W⋀U]

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[1]⋀S0[1]⋀N1[2]]

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 (see section 2.4) have been defined so as to take into account the missing identity matrix.

2.3  Auxiliary functions / global variables

(i) FormDegree, simpRules, varList, zeroQ

•	The function FormDegree[x_] returns the (differential-form) degree of the expression x.  It assumes that all terms in a sum and all components in an explicit matrix have the same degree.

FormDegree[a1⋀b2]

FormDegree[S⋀U]

FormDegree[d[S]⋀U]

FormDegree[R3]

•	The global variable simpRules is a list of "simplification rules" that the user can define. This list is used (repeatedly) by reWrite to simplify expressions.  The default value of  simpRules  is: simpRules = {x_.  Cos[z_]^2 + x_.  Sin[z_]^2 :> x}. It's effect can be seen in the following example:

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

reWrite[%]

•	Another global variable, useful in simplifying complicated expressions, is varList. It is a list of variables (or Patterns) that the user can define; they are used by reWrite  to group terms in the coefficients of identical Wedge[__] or  d[_] terms, so as to further simplify the result. The default value of  varList  is the empty list: varList = {}. An illustration of its use is given in the following example:

Example of 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]

which is a more meaningful expression. Verify that it is a closed form:

reWrite[d[%]]

varList = {};

Note: The global variables simpRules and varList are assumed to be lists. Thus, to clear them, one must not use Clear[simpRules, varList].  One must instead set them equal to the empty list.

•	The function zeroQ[x_] tests if the expression x  is identically zero (zero matrix) and returns True  or False.  It does not perform any simplifications:

zeroQ/@{0, a, Cos[x]^2 + Sin[x]^2 - 1, {{0, 0}, {0, 0}}, 
        a^2 + 2  a b + b^2 - (a + b)^2}

It should be used after simplification with reWrite:

zeroQ/@
      reWrite[{0, a, Cos[x]^2 + Sin[x]^2 - 1, {{0, 0}, {0, 0}}, 
          (a + b)^2 - a^2 - 2  a b - b^2}]

(ii) Bar

The function Bar[x_] returns the complex conjugate of the expression x, assuming all symbols are complex.

f2

Bar[%]

If some symbols are real, the user must either set Bar[ai]=ai, or use a corresponding replacement rule, for each real symbol ai.

% /. {Bar[x] -> x, Bar[y] -> y, Bar[Cos] -> Cos, 
        Bar[Sin] -> Sin}

If most symbols are real, as in this example, one can define Bar[x_Symbol]:=x,  and still set Bar[ci]=ciB; Bar[ciB]=ci; for pairs of complex-conjugate symbols ci, ciB (immediate assignments are evaluated before delayed ones!).

r1 = {z -> x + I y}

r2 = Join[r1, Bar[r1]]

\(Bar[x_Symbol] := x; Bar[z] = zB; Bar[zB] = z; 
r2\)

Bar[%]

Bar[x_Symbol] =.

(iii) ScalFormQ, BasicScalFormQ, MatQ, BasicMatQ

•  The functions 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. Note that the complex conjugate of a basic form and the trace of a p(>0)-form matrix are taken to be basic forms (the symbols below have been defined as symbolic forms or matrices in section 2.2).

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

ScalFormQ/@testExpr

BasicScalFormQ/@testExpr

•  The functions 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], a1, S, d[S], N1[2] + N1t[2], S⋀U, d[x]⋀U, 
        Bar[N1[2]]};

MatQ/@testExpr

BasicMatQ/@testExpr

These functions are useful in setting up conditions for applying extra Wedge definitions or rules.

(iv) NoDif

See Example 3.5

(v) indepTerms, indepCoefs

These functions are useful in reducing large lists of expressions.

•	The function indepTerms[x_]  returns a list of the "independent" elements of the list x. Independent here means not equal to  ±  a numerical multiple of another element. Note that indepTerms does not perform any reductions:

indepTerms[{2  a + 2  b, a + b}]

indepTerms[{2 (a + b), a + b}]

nor does it test for linear dependence:

indepTerms[{a - b, b - c, c - a}]

Plus@@%

•	The function indepCoefs[x_]  returns a list of the "independent" coefficients of all symbolic differential forms, symbolic matrices or wedge-products in the elements of the list x. Independent here means not equal to  ±  a numerical multiple of another coefficient.

The symbols q[i] have been defined earlier as symbolic 1-forms.

indepCoefs[a[1] q[1] + a[2] q[2] + a[3] q[3]]

indepCoefs works with symbolic matrices also (S has been defined earlier as a 0-form symmetric matrix).

indepCoefs[
      a[0] + a[1] S + a[2] S⋀S + a[3] S⋀S⋀S]

indepCoefs[
      (a[2] + a[3]) q[2]⋀q[3] + 
        (a[3] + a[1]) q[3]⋀q[1] + 
        (a[1] - a[2]) q[1]⋀q[2]]

Like indepTerms, indepCoefs does not test for linear dependence:

Plus@@%

indepCoefs has no effect on expressions that are not differential forms or matrices:

indepCoefs[a[1] b[1] + a[2] b[2] + a[3] b[3]]

indepCoefs can be applied to lists:

sym33mat = 
      Table[(i + j) (a[i] - a[j]) q[i]⋀q[j], {i, 3}, {j, 3}]

indepCoefs[%]

Note that indepTerms does not recognize that (-a[1]+a[2])=-(a[1]-a[2]), as it does not carry out any simplifications:

indepTerms[sym33mat]

indepTerms[reWrite[sym33mat]]

(vi) GenCoef

The function  GenCoef[x_,y_] ("Generalized Coefficient") is very similar to the Mathematica function Coefficient[x_, y_]. Unlike Coefficient, however, GenCoef returns a series when x is a series and, when both x and y are Lists, returns a nested List:

seriesExpr = Series[Sum[(a + b)^k*x^k, {k, 0, 3}], {x, 0, 3}]

Coefficient[seriesExpr, b]

Coefficient[seriesExpr, b, 2]

GenCoef[seriesExpr, b]

GenCoef[seriesExpr, b, 2]

Coefficient[{a[1] x + a[2] y, b[1] x + b[2] y}, {x, y}]

GenCoef[{a[1] x + a[2] y, b[1] x + b[2] y}, {x, y}]

Thus GenCoef[LinearEquations_List,unknowns_List] gives directly the matrix of coefficients of the unknowns in the given equations.

(vii) FormCoef

The function  FormCoef[x_,y_] ("Form Coefficient") gives the "left-coefficient" of the differential form y in the differential form expression x, i.e., writes x as w⋀y + terms not containing y and returns w. If y is of odd FormDegree, then w is defined up to a multiple of y. The differential form y must either be a single symbolic form or a Wedge product of such forms. Examples (the symbols e[i] are defined as symbolic 1-forms):

threeForm = 
      d[f]⋀e[1] ⋀e[3] + 
        (f + g) e[1]⋀ e[2]⋀e[3] + 
        d[h]⋀e[1] ⋀e[2] + h e[1]⋀d[e[2]]

FormCoef[threeForm, e[1]]

FormCoef[threeForm, e[2]]

FormCoef[threeForm, d[e[2]]]

FormCoef[threeForm, d[f]]

FormCoef[threeForm, e[1]⋀e[2]]

FormCoef[threeForm, d[h]⋀e[2]]

FormCoef[threeForm, e[1]⋀e[2]⋀e[3]]

FormCoef[threeForm, d[f]⋀e[1] ⋀e[3]]

FormCoef[threeForm, f e[1]⋀e[2]]

FormCoef gives here the wrong result because f e[1]⋀e[2] is neither a single symbolic form nor a Wedge product of such forms.

The arguments of FormCoef can contain symbolic forms of any degree: 
	(the symbols {a2, b2, c2} have been defined earlier as symbolic 2-forms)

fourForm = 
      e[1]⋀ e[2] ⋀a2 + e[2]⋀e[3]⋀ b2 + 
        e[3] ⋀e[1]⋀c2 + a2⋀b2 + b2⋀c2 + 
        a2⋀c2

FormCoef[fourForm, a2]

FormCoef[fourForm, e[1]]

FormCoef[fourForm, e[1]⋀a2]

If x is a list of differential forms, FormCoef returns the corresponding list of coefficients:

FormCoef[{e[1]⋀e[2], e[2]⋀e[3], e[3]⋀e[1]}, e[1]]

2.4  Extra matrix functions / global variables

Even though symbolic matrices can have unequal dimensions, some operations defined in this section (tests for symmetry, trace, etc.) can only be applied to square matrices.

(i) Transpose, SymQ, AntisymQ

•	As mentioned in section 2.2 (ii), the built-in function Transpose  has been extended to handle properly symbolic matrix expressions containing the symbols introduced with DeclareMatrixForms. (The symbols S and U, A1[_] have been defined above as a symmetric 0-form and antisymmetric 1-form matrices.)

Transpose[S⋀U+U⋀S]

d[S⋀U+U⋀S]+A1[1]⋀(S⋀U+U⋀S)+(S⋀U+U⋀S)⋀A1[1]

% + Transpose[%]

Using Transpose, the following two functions have been defined:

•	The function SymQ[x_] returns True  if  Transpose[x]=x and False otherwise.
The function 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[%]

If x is a scalar, it is assumed to multiply the identity matrix, so SymQ/AntisymQ  give True /False

SymQ[a]

AntisymQ[a]

(ii) matDimension, traceOption, tr

•	The global variable matDimension is a positive integer -- the matrix dimension (assumed square) -- to be assigned by the user. All symbolic matrices are assumed to have the same dimension. Together with a second global variable, traceOption, it is used in symbolic trace calculations (see below).

•	The function tr[x_] returns the trace of the (symbolic or explicit, square) matrix expression x. If  x is a scalar, it is assumed to multiply the identity matrix, and thus tr[x] 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;(** No rules applied, other than linearity and tr[At] = tr[A] **)tr/@testExpr

traceOption = 1;(** Cyclic rule also applied **)tr/@testExpr

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

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 (alphabetical) 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}

(iii) toComponents

The function toComponents[x_,rule_List] takes a symbolic matrix expression (x) and a list of rules (rule)  and returns the explicit form (list of lists) of the symbolic matrix expression, introducing the identity matrix where needed. The explicit form of all basic symbolic matrices in x must be given as a set of replacement rules in the list rule 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  rule gives the components of U  (which has been declared an antisymmetric 1-form matrix) in terms of the undefined scalar 1-forms u1, u2, u3, one must evaluate DeclareForms[1,u1,u2,u3] before using toComponents.

matDimension = 3;
matCompRules={A -> Array[aa, {3, 3}],At -> Transpose[Array[aa, {3, 3}]],S -> Array[ss, {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):

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

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

DeclareForms[1, uu[_]]

matCompRules

MatrixForm[toComponents[S⋀U,matCompRules]]

tr[%]

FormDegree[%%]

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

MatrixForm[toComponents[A  +  lamda, matCompRules]]

tr[%]

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

tr[%]

MatrixForm[toComponents[%%, matCompRules]]

SymQ[%]

tr[%%]

Collect[%, lamda, Factor]

matDimension =.

(iv) CHrule, detTOtr

•	The function CHrule[n_Integer, mat_:$] returns a rule for implementing the Cayley-Hamilton theorem for the  0-form symbolic matrix mat  with  matDimension=n. The second argument is optional; if it is missing CHrule returns a delayed rule that can be applied to any n x n  0-form matrix, using the symbol $ which is defined in EDC as a 0-form symbolic matrix. 
Note: After calling CHrule[n], matDimension is set to n.

CHrule[3, A]

matDimension

CH3 = CHrule[3]

The rule applies only when $ is a 0-form matrix expression. It is applied for products of matrices also:

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

CH4 = CHrule[4]

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

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

For antisymmetric matrices the trace function (tr) gives 0 for odd powers, simplifying the expression:

DeclareMatrixForms[{0, H, -H}]

CHrule[4, H]

The same result is obtained using the general rule CH4:

%[[1, 1]] /. CH4

Warning : The rules returned by CHrule 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 //. CHrule[3, A], 2]

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

•	The function detTOtr[n_Integer, mat_:$] uses CHrule to express the determinant of the   0-form symbolic matrix  mat  with  matDimension=n, as an n-th order polynomial in the traces of  mat. When the second argument is omitted, the symbol $ is used for the matrix. Examples:

detTOtr[2]

detTOtr[3, A]

detTOtr[4, B]

Using explicit matrices, one can check that detTOtr gives the correct result:

detTOtrTEST[n_] := 
      (nxnMat = Array[mij, {n, n}]; 
        Factor[detTOtr[n] /. $ -> nxnMat] === Factor[Det[nxnMat]])

detTOtrTEST/@{2, 3, 4, 5}

(v) subMatTK, subMatDR

The functions subMatTK[x_, rr_, cc_] and subMatDR[x_, rr_, cc_] return the sub-matrices that are obtained by Taking (TK) or Dropping (DR) the rows indicated by rr and the columns indicated by cc of the matrix x. The parameters  rr and cc have the same form as in the corresponding Mathematica commands (Take, Drop): either a list of two integers {n1, n2} indicating the range {n1...n2} of rows (columns), or a single integer n which is interpreted as {1, n}. For example, let mat67 be a 6x7 matrix:

mat67 = Array[mm, {6, 7}]

The following matrix consists of the entries in rows (2,3,4) and columns (6,7) of mat67:

subMatTK[mat67, {2, 4}, {6, 7}]

The following matrix consists of the entries in rows (1,5,6) and columns (1,2,3,4,5) of mat67:

subMatDR[mat67, {2, 4}, {6, 7}]

Similarly, taking or dropping the first 3 rows and columns (5,6) of mat67:

subMatTK[mat67, 3, {5, 6}]

subMatDR[mat67, 3, {5, 6}]

(vi) DirectSum, FlatOuter

Let us define some explicit matrices:

A22 = Array[a, {2, 2}]

B33 = Array[b, {3, 3}]

P32 = Array[p, {3, 2}]

Q23 = Array[q, {2, 3}]

MatrixForm/@{A22, B33, P32, Q23}

•	The function DirectSum[x__] takes two or more explicit matrices of dimensions (Subscript[r, i]) x (Subscript[c, i]) and constructs their direct sum, a matrix with dimensions  Σ (Subscript[r, i]) x Σ (Subscript[c, i]), having the given matrices along the main diagonal and 0 everywhere else:

MatrixForm[DirectSum[A22, B33, P32, Q23]]

Dimensions[%]

Plus@@(Dimensions/@{A22, B33, P32, Q23})

•	The function FlatOuter[x__] takes two or more explicit matrices of dimensions r_i x c_i and constructs their outer product (KroneckerProduct), a matrix with dimensions  Π (Subscript[r, i]) x Π (Subscript[c, i])  (KroneckerProduct is a built-in function in Mathematica 6 or later):

MatrixForm[FlatOuter[A22, P32]]

MatrixForm[FlatOuter[P32, A22]]

MatrixForm[FlatOuter[P32, A22, Q23]]

Dimensions[%]

Times@@(Dimensions/@{P32, A22, Q23})

MatrixForm[FlatOuter[P32, Q23, A22]]

FlatOuter[P32, Q23, A22] === FlatOuter[FlatOuter[P32, Q23], A22]

FlatOuter[P32, Q23, A22] === FlatOuter[P32, FlatOuter[Q23, A22]]

If[$VersionNumber > 5.5, 
      FlatOuter[P32, A22, Q23] === KroneckerProduct[P32, A22, Q23]]

(vii) Contract, multiDot

These functions perform contractions on nested lists. It is convenient to think of an nth-level nested list as an nth-rank tensor. Contraction then produces lower rank tensors. For use in the examples we define the following rank-3 and rank-4 tensors in three dimensions:

\(T3 = Array[t3, {3, 3, 3}]; 
T4 = Array[t4, {3, 3, 3, 3}]; \)

•	The function Contract[x_,{i1,j1},{i2,j2},...] takes as arguments an nth-rank tensor x, and one or more lists of pairs of integers, indicating the positions of the indices to be contracted. It returns a tensor of rank n-2k, where k is the number of lists:

Contract[T3, {1, 2}]

Contract[T3, {1, 3}]

Contract[T4, {1, 2}]

Contract[T4, {1, 2}, {3, 4}]

Contract[T4, {1, 4}, {2, 3}]

•	The function multiDot[x_,y_,{i1,j1},{i2,j2},...] is similar to Contract, and generalizes the built-in function Dot. It takes as arguments two tensors -- x, y -- and one or more lists of pairs of integers. For each such pair {i, j}, it contracts the ith index of the first tensor (x) with the jth index of the second tensor (y), returning a tensor of rank m+n-2k, where m, n are the ranks of  x, y, respectively, and k is the number of index pairs. (When the components of x and y are differential forms, multiDot multiplies them using Wedge).

multiDot[T3, T3, {1, 1}]

Dimensions[%]

multiDot[T3, T3, {1, 1}, {2, 2}]

Dimensions[%]

multiDot[T3, T3, {1, 2}, {2, 1}]

multiDot[T3, T3, {1, 1}, {2, 2}, {3, 3}]

If we declare the t3[__] to be 1-forms, this contraction vanishes, as multiDot uses Wedge:

DeclareForms[1, t3[_]]

multiDot[T3, T3, {1, 1}, {2, 2}, {3, 3}]

multiDot[T3, T3, {1, 2}, {2, 1}]

multiDot[T3, T4, {1, 1}, {2, 2}, {3, 3}]

Dimensions[%]

multiDot[T4, T4, {1, 1}, {2, 2}, {3, 3}]

Dimensions[%]

Note1: In both Contract and multiDot "contraction" is simple summation; the tensor indices to be contracted must have opposite variance (one co- and one contra- variant).

Note2: One can use Contract and multiDot to contract nested lists with unequal index dimensions (provided, of course, that the indices to be contracted have the same dimension); thus, if A has dimensions {2, 3, 3, 4} and B {3, 2, 4}, Contract[A, {2,3}] will return a tensor with dimensions {2, 4} and multiDot[A, B, {1,2}, {3,1}] will return a tensor with dimensions {3, 4, 4}.

2.5  interiorProduct, LieDcartan

Two new functions introduced in version 3.8.5 are interiorProduct[vector_,form_], giving the interior product (contraction) between a vector and a p-form, resulting in a (p-1)-form, and LieDcartan[vector_,form_], giving the Lie derivative of the differential form expression form w.r.t. the vector. LieDcartan is defined in terms of interiorProduct using Cartan's identity. interiorProduct (and hence LieDcartan) only applies to p-forms that are expressed as wedge products of certain basic 1-forms. It has no effect on symbolic p-forms.

The vector in the first argument of interiorProduct[vector_,form_], must be given as a list of its contravariant components in a frame dual to the basic 1-forms used: if the coframe basis is, say, {e[1], e[2], ..., e[n]}, then the list {(V^1),(V^2),...,(V^n)} will denote the vector  (V^1)(Subscript[X, 1])+ (V^2)(Subscript[X, 2])+... +(V^n)(Subscript[X, n]), where the frame basis vectors (Subscript[X, k]) satisfy (Subscript[X, a])(e^b)=((Subscript[δ, a])^b) ((e^b)= e[b]). The symbols e[_] are declared as 1-forms in the matrixEDC code.

The coframe basis used to express all p-forms must be either a numbered set like {σ[1], σ[2], ..., σ[n]} or {ω[1], ω[2], ..., ω[n]} (the symbols σ[_], ω[_] must first be declared as 1-forms), or the differentials of the coordinates {d[x], d[y],...,d[w]}. However, when the coordinate differentials are used interiorProduct requires a third argument -- the ordered list of coordinates -- so as to know which vector component corresponds to which coordinate. For example, 
interiorProduct[{1,0,0,0}, d[x]⋀d[t],{x,y,z,t}] gives d[t], but 
interiorProduct[{1,0,0,0}, d[x]⋀d[t],{t,x,y,z}] gives -d[x].

The following examples illustrate the use of interiorProduct and LieDcartan:

vector = {v[1], v[2], v[3], v[4]}

test1form = X e[1] + Y e[2] + Z e[3] + T e[4]

test2form = 
      P  e[1]⋀e[4] + Q e[2]⋀e[4] + R  e[3]⋀e[4] + 
        E e[2]⋀e[3] + F  e[3]⋀e[1] + G e[1]⋀e[2]

test3form = 
      K  e[1]⋀e[2]⋀e[3] + L e[2]⋀e[3]⋀e[4]  + 
        M  e[1]⋀e[3]⋀e[4] + 
        N e[1]⋀e[2]⋀e[4]

interiorProduct[vector, test1form]

interiorProduct[vector, test2form]

interiorProduct[vector, test3form]

The first argument of interiorProduct must be a list:

interiorProduct[test1form, test3form]

interiorProduct has no effect on symbolic forms:

DeclareForms[{1, form1},  {2, form2}]

interiorProduct[vector, form2]

•	interiorProduct can also be applied to lists of forms, giving the corresponding list of contractions with the given vector:

interiorProduct[
      vector, {e[1]⋀e[2]⋀e[3], e[2]⋀e[3]⋀e[4], 
        e[3]⋀e[4]⋀e[1], e[4]⋀e[1]⋀e[2]}]

form2array = Array[Wedge[e[#1], e[#2]]&, {4, 4}]

interiorProduct[vector, form2array]

•	One can use any numbered set of 1-forms as basis 1-forms:

DeclareForms[1, σ[_], ω[_]]

test2form /. e -> σ

interiorProduct[vector, %]

test3form /. e -> ω

interiorProduct[vector, %]

•	However, when the coordinate differential basis is used, a third argument is required to define the order of the vector components:

interiorProduct[{1, 0, 0, 0},  d[x]⋀d[t]]

interiorProduct[{1, 0, 0, 0},  d[x]⋀d[t], {x, y, z, t}]

interiorProduct[{1, 0, 0, 0},  d[x]⋀d[t], {t, x, y, z}]

interiorProduct[{1, 0, 0, 0},  d[x]⋀d[t], {z, x, y, t}]

•	The Lie derivative of forms, LieDcartan, has been defined using interiorProduct and Cartan's identity. But as it involves exterior derivatives, which treat all symbols as variables, it cannot be completely evaluated unless one gives the exterior derivatives of the basis 1-forms, and of all other symbols used, in terms of the basis 1-forms. For example:

d[e[4]] = 0;

d[e[i_]] := 
      1/2  Sum[Signature[{i, j, k}] Wedge[e[j], e[k]], {j, 3}, {k, 3}]

LieDcartan[vector, test2form]

To complete the evaluation, d[v[k]], d[P], d[Q], etc., must be given in terms of the basis e[a]. For example, consistent with our choice of d[e[i]], we can define:

d[v[4]] = e[4]; d[v[k_]] = 0;

d[x_Symbol] := 0 /; FormDegree[x] === 0

LieDcartan[vector, test2form]

As with interiorProduct, LieDcartan requires a third argument in a coordinate frame:

d[x_Symbol] =. ; d[v[4]] =. ;

LieDcartan[vector, 
      K[x, y, z] d[x]⋀d[t] + L[x, y, z] d[y]⋀d[t] + 
        M[x, y, z] d[z]⋀d[t]]

LieDcartan[vector, 
      K[x, y, z] d[x]⋀d[t] + L[x, y, z] d[y]⋀d[t] + 
        M[x, y, z] d[z]⋀d[t], {x, y, z, t}]

Note that, because of its definition in terms of interiorProduct, LieDcartan acts on lists as a collection of scalar p-forms and not as the components of a tensor. In the combined EDCRGTCcode a second Lie derivative, LieD, is defined that treats lists as the components of a tensor.

The operators d and LieDcartan commute, even though the calculations cannot be completed, due to d[d[anything]]=0:

reWrite[
      d[LieDcartan[vector, test2form]] - LieDcartan[vector, d[test2form]]]

reWrite[
      d[LieDcartan[vector, test3form]] - LieDcartan[vector, d[test3form]]]

For the same reason, even though LieDcartan does not evaluate on symbolic forms, its commutator with d vanishes:

LieDcartan[vector, form2]

reWrite[d[LieDcartan[vector, form2]] - LieDcartan[vector, d[form2]]]

d[v[k_]] =. ; d[e[i_]] =. ;

2.6  Input Palette

The  notebook EDCpalette.nb is an input palette containing most of the new functions/variables, so they can be entered easily without misspellings.

3  Examples

All symbolic calculations with EDC begin with (i) defining symbols for the basic, undefined forms and matrices one intends to use and (ii) giving definitions for the action of Wedge, d, tr, etc. on these symbols. The following examples illustrate these steps. More examples in the notebook  matrixEDCexamples.nb.

3.1  Symbolic Identity Matrix

Allthough EDC uses the number 1 in symbolic matrix expressions for the identity matrix, one can give appropriate definitions to make any other symbol behave as a symbolic identity matrix. For example:

DeclareMatrixForms[{0, IdMat, IdMat}]

tr[IdMat] = matDimension;

d[IdMat] = 0; Bar[IdMat] = IdMat;

Wedge[IdMat, IdMat] = IdMat;

Wedge[IdMat, x_?MatQ] := x; Wedge[x_?MatQ, IdMat] := x;

These rules are used below:

\(Alam = A - lamda*IdMat; 
Wedge[Alam, Alam, Alam]\)

d[%]

3.2  Polynomial Matrix Algebra

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

DeclareMatrixForms[{0, Z[_], Z[_]}]

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

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

CH3 = CHrule[3, Z[1]]

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

For simplicity (and to avoid infinite loops!) set

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

CH3rhs

In this notation Z[0] is the identity matrix. One can either set Z[0]=1,  or set Z[0]=IdMat (defined in Example 3.1). For this example however, as only products with Z[i] are involved, it is enough to set

tr[Z[0]] = 3;

and use  Z[0]  to rewrite CH3rhs:

CH3rhs = reWrite[CH3rhs⋀Z[0]]

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

Now define rules for simplifying Z[i] matrices

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

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

Z[10] // Timing

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

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

reWrite[matPoly⋀matPoly⋀matPoly] // Timing

Finding the inverse of  matPoly

Take a linear combination of the powers of the given matrix, and determine the coefficients by requiring that it is the inverse of the matrix polynomial matPoly:

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

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

eqs = {zer0 /. Z[0] -> 1, zer0 /. Z[1] -> 1, zer0 /. Z[2] -> 1} /. 
        Z[_] -> 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.3  SU(2) Trace Calculations

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

Introduce notation and define properties of the 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

Clear[matDimension, A, B]; 
    On[General::"spell", General::"spell1"];

3.4  Non-commutative Algebra: Applying Commutators

A common type of calculation involves evaluation of polynomial expressions in non-commuting variables. For such calculations one introduces symbolic matrices for these variables and gives appropriate Wedge rules. We will illustrate this with a pair of canonically conjugate quantum mechanical variables P, Q, satisfying [Q, P]=I h. As in Example 3.2 we will use the notation P[n] for P^n.

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

DeclareMatrixForms[{0, P[_], Pt[_]}, {0, Q[_], Qt[_]}]

P[0] = 1; Q[0] = 1;

Wedge[P[m_], P[n_]] := P[m + n]; Wedge[Q[m_], Q[n_]] := Q[m + n];

Wedge[Q[m_], P[n_]] := 
      Wedge[Q[m - 1], Wedge[P[1], Q[1]] + I h, P[n - 1]]

The last definition implements the commutator; it will be applied repeatedly to put P^n before Q^m. 
(To put Q's before P's use:  Wedge[P[m_],Q[n_]]:=Wedge[P[m-1],Wedge[Q[1],P[1]]-I h,Q[n-1]]).

Wedge[P[2], Q[2]]

Wedge[Q[2], P[2]]

reWrite[%]

We will now define 3 polynomials in P,Q with arbitrary coefficients:

poly1 = Sum[a[m, n] P[m]⋀Q[n], {m, 0, 3}, {n, 0, 3}]

poly2 = Sum[b[m, n] P[m]⋀Q[n], {m, 0, 3}, {n, 0, 3}]

poly3 = Sum[c[m, n] P[m]⋀Q[n], {m, 0, 3}, {n, 0, 3}]

One can then compute powers or commutators of these polynomials (these calculations take long!):

reWrite[poly1⋀poly1]

comm[x_, y_] := reWrite[x⋀y - y⋀x]

comm[poly1, poly2]

Collect[%, h, reWrite]

JacobiID[x_, y_, z_] := 
      reWrite[comm[x, comm[y, z]] + comm[y, comm[z, x]] + comm[z, comm[x, y]]]

JacobiID[poly1, poly2, poly3] // Timing

One can speed up the calculation by giving extra rules for commuting higher powers of P, Q:

Q2P2 = reWrite[Wedge[Q[2], P[2]]]

Q3P3 = reWrite[Wedge[Q[3], P[3]]]

Q2P3 = reWrite[Wedge[Q[2], P[3]]]; Q3P2 = reWrite[Wedge[Q[3], P[2]]];

Wedge[Q[m_], P[n_]] := 
      Wedge[Q[m - 3], Q3P3, P[n - 3]] /; m > 2 && n > 2

Wedge[Q[m_], P[n_]] := 
      Wedge[Q[m - 3], Q3P2, P[n - 2]] /; m > 2 && n > 1

Wedge[Q[m_], P[n_]] := 
      Wedge[Q[m - 2], Q2P3, P[n - 3]] /; m > 1 && n > 2

Wedge[Q[m_], P[n_]] := 
      Wedge[Q[m - 2], Q2P2, P[n - 2]] /; m > 1 && n > 1

JacobiID[poly1, poly2, poly3] // Timing

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

3.5  Introducing Custom Notation: Symbolic Inverse

Sometimes one may want to use a shorthand notation of the form K[A] for a symbolic operation with matrices (or other variables) that should not be interpreted as a function by the exterior differential operator d; i.e., one wants to avoid the following:

d[{Inverse[S], K[A, B]}]

•	The command  NoDif[heads__]  ("No Differentiation") takes as arguments the "Heads" of the expressions of the form  K[A]  that must not be differentiated and blocks the action d:

NoDif[Inverse, K]

Now d  has no effect:

d[{Inverse[S], K[A, B]}]

Next, one must enable EDC to decide if these expressions represent differential forms or matrices. This is achieved by:
	(1) Setting the value of FormDegree[..] if > 0 and
	(2) Setting the values of BasicMatQ[..] / BasicScalFormQ[..] if True. (If BasicMatQ[x] gives True, then MatQ[x] gives True also.)

For example,

\(BasicMatQ[Inverse[x_]] := BasicMatQ[x]; 

    FormDegree[Inverse[x_]] := FormDegree[x]; 

    FormDegree[K[x_, y_]] := FormDegree[x]; \)

MatQ[Inverse[A]]

FormDegree[K[U, A]]

FormDegree[Inverse[A]]

One can then give definitions for the action of d and Wedge on these expressions as well as other definitions. The following definitions enable the built-in function Inverse to give the symbolic inverse of symbolic matrices (on matrix products only! These definitions have not been included in EDC as they cannot be extended to sums of symbolic matrices):

\(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; \)

Unprotect[Inverse]; Inverse[Inverse[x_]] = 1; 
    Inverse[a_*x_] := 1/a*Inverse[x] /; FormDegree[a] === 0 && ! MatQ[a]; 
    Inverse[Wedge[x__]] := 
      Wedge@@(Reverse[Inverse/@{x}]) /; Union[BasicMatQ/@{x}] === {True}; 
    Protect[Inverse];

Inverse[A⋀B⋀C]

Wedge[A⋀B⋀C, %]

Wedge[C, %%, A]

S⋀d[Inverse[S]]

d[Inverse[S]]⋀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;!

3.6  Symmetries of general 2nd order ODE

We first reevaluate the initialization code matrixEDC.m to clear previous symbol definitions, as we will use some of the same symbols with a different meaning.

<< matrixEDC.m

We then define the function FuncRepRules (which is defined in RGTC) that we will need for function replacement:

FuncRepRules[z_[y__], g_, ord_:  2] := 
      FuncRepRules[z[y], g, ord] = 
        Block[{tmp$ = {z[y] -> g}, indL = indexList[Length[{y}], ord]}, 
          Do[tmp$ = {tmp$, 
                Derivative[Sequence@@indL[[i]]][z][y] -> 
                  Expand[FaSi[
                      D[g, Sequence@@Transpose[{{y}, indL[[i]]}]]]]}, {i, 
              Length[indL]}]; Flatten[tmp$]];

\(\(indexList[varNo_, ord_] := 
      Block[{tmp$, indL}, 
        tmp$ = NestList[RotateRight, Join[{1}, Table[0, {varNo - 1}]], 
            varNo - 1]; 
        If[ord > 1, 
          Do[indL = tmp$; 
            Do[tmp$ = Join[tmp$, Map[Plus[indL[\([i]\)], #]&, indL]], {i, 
                varNo}], {ord - 1}]]; Union[tmp$]]
\)\)

The general 2nd order ODE   u''[x]-f[x,u,u'[x]]==0  is equivalent to the Exterior Differential System (EDS) generated by the pair of 1-forms:

EDS = {d[u] - p d[x], d[p] - f[x, u, p] d[x]}

The pair of 1-forms EDS generate a closed differential ideal:

mat22 = Array[m, {2, 2}]

reWrite[d[EDS] + mat22.EDS⋀d[x]]

reWrite[d[EDS]+{{0, 1},{Derivative[0, 1, 0][f][x, u, p],Derivative[0, 0, 1][f][x, u, p]}}.EDS⋀d[x]]

A symmetry of the system is generated by a vector field with respect to which the Lie derivative of the EDS vanishes (modulo the EDS). This implies that:

vectorField = {X[x, u, p],  U[x, u, p], P[x, u, p]};

matEQS = 
      reWrite[LieDcartan[vectorField, EDS, {x, u, p}] - mat22.EDS]

GenCoef[matEQS, {d[p], d[u]}]

matSOL = Solve[% == 0, Flatten[mat22]][[1]]

EDSeqs = reWrite[matEQS /. matSOL]

A solution is obtained if the vectorField is assumed to be proportional to the total derivative w.r.t. x, i.e., U[x,u,p]=p X[x,u,p] and P[x,u,p]=f[x,u,p] X[x,u,p]:

reWrite[
      EDSeqs /. FuncRepRules[U[x, u, p], p X[x, u, p]] /. 
        FuncRepRules[P[x, u, p], f[x, u, p] X[x, u, p]]]

The function P[x,u,p] does not appear differentiated in the first of EDSeqs, so it can be expressed in terms of the other functions:

Psol = Solve[EDSeqs[[1]] == 0, P[x, u, p]][[1]]

lastEQ = 
      indepTerms[
          reWrite[EDSeqs /. FuncRepRules[P[x, u, p], Psol[[1, 2]]]] /. 
            d[x] -> 1][[1]]

This last equation depends on 3 arbitrarty functions (X[x,u,p], U[x,u,p], f[x,u,p]). A particular solution, depending on an arbitrary function of one variable (g[p]), is X[x,u,p]=1, U[x,u,p]=1, and f[x,u,p]=(x-u)g[p] :

Factor[
      lastEQ /. FuncRepRules[X[x, u, p], 1] /. 
          FuncRepRules[U[x, u, p], 1] /. 
        FuncRepRules[f[x, u, p], (x - u) g[p]]]

It can be verified that, for these choices of X[x,u,p], U[x,u,p], and f[x,u,p], the function P[x,u,p] vanishes. Clearly there are many more solutions.

EDC is also used in the Riemannian Geometry & Tensor Calculus package RGTC (Sections 5 and 6). There is also an extension of EDC, called superEDC, that allows calculations with graded (Grassmann) variables.

4  Usage Tips

4.1  Explicit vector calculations

Wedge  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 = Array[a, {3, 3}]

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⋀%]

4.2  Using lists of symbolic matrices

Occasionally, one may want to do calculations with lists of symbolic matrices. For example, consider a 2x2 matrix (or vector) in block-form:

DeclareMatrixForms[{0, M[_], Mt[_]}]

blockMat = {{M[11], M[12]}, {M[21], M[22]}}

blockVec = {M[1], M[2]}

Wedge does not multiply such lists with a single symbolic matrix, but it does multiply them together:

Wedge[M[1], blockVec]

Wedge[blockMat, M[1]]

Wedge[blockMat, blockVec]

Wedge[blockMat, blockMat]

Wedge[blockVec, blockMat]

The reason for this strange behavior is compatibility of dimensions: if the M[i] are n x n matrices, then a single matrix M[1] cannot multiply a 2n x 2n matrix (blockMat) or a 2n vector (blockVec).  But 2n-dim matrices and vectors can be multiplied together.

However, if one's symbolic matrices have a different meaning and components will never be used, one can force Wedge to perform such multiplications. This is achieved by evaluating the following Wedge definitions:

Wedge[x_, y_List] := Map[Wedge[x, #]&, y]

Wedge[x_List, y_] := Map[Wedge[#, y]&, x]

Now:

Wedge[M[1], blockVec]

Wedge[blockMat, M[1]]

4.3  Calculations with Series Expansions

When working with series expansions containing explicit differential forms, 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 the expansion variable (and appropriately define its exterior derivative). For example, the following definitions introduce the symbol dx to be used in place of d[x]:

\(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

d[x] =.