We begin by turning off Mathematica's spelling errors, since we will be using many variables which are spelled similarly:
Off[ General::spell, General::spell1];We now use the SetDelayed operator, previously used to define Mathematica functions, to define the algebraic properties of a generalized derivative operator without defining the actual operator itself:
Clear[ der];We began by using the Mathematica Clear function to remove any previous definitions for the symbol "der". This is important because when defining a function in this manner, the order in which the properties are specified controls how the function is evaluated: the "rules" are applied in the order in which they were defined. We then specified that the derivative of a numeric quantity is zero, and explained to Mathematica how the derivative operator distributes over addition and multiplication.der[ expr_?NumericQ, x_] := 0;
der[ expr_Plus, x_] := Sum[ der[ expr[[iind]], x], {iind, Length[expr]}];
der[ expr_Times, x_] := Sum[
Product[ expr[[jind]], {jind, iind - 1}] * der[ expr[[iind]], x] *Product[ expr[[jind]], {jind, iind + 1, Length[expr]}], {iind, Length[expr]}];
We repeat the process for the variational operator (usually denoted δ):
Clear[ var];We further specify that the variation of a derivative is the derivative of the variation:var[ expr_?NumericQ] := 0;
var[ expr_Plus] := Sum[ var[ expr[[iind]]], {iind, Length[expr]}];
var[ expr_Times] := Sum[
Product[ expr[[jind]], {jind, iind - 1}] * var[ expr[[iind]]] *Product[ expr[[jind]], {jind, iind + 1, Length[expr]}], {iind, Length[expr]}];
var[ der[ expr_, x_]] := der[ var[expr], x];Note that we do not define the properties of the covariant derivative; we will not need them, and in fact it will prove useful below to express the covariant derivative of a product without performing the distribution.var[ covder[ expr_, x_]] := covder[ var[expr], x];
Many of the arguments to the functions we will be working with will be indices. It will therefore be
necessary to generate unique index variables for summations. While Mathematica's Unique function
will do this, the resultant expressions are less than clear for anyone used to the usual notation.
We therefore define a function ind which has no arguments and whose value is a symbol constructed
from two or more identical letters (keeping with our practice of not using single character variables
in functions):
If[ indi > 26, index = Mod[indi - 1, 26] + 1, index = indi];
num = Quotient[indi + 25, 26] + 1;
str = "";
Do[ str = str <> abet[[index]], {ii, num}];
Return[ Symbol[str]];)];
We will be analyzing expressions which contain a number of both free and contracted indices. In many
cases, simplification of these expressions will be very hard for Mathematica because it does not know
to rename indices in order to facilitate factoring. We therefore define a function iclean, which
will attempt to clean up contracted indices and do preliminary simplification. Note that there
is no easy way to do this in a foolproof manner; some manual simplification will inevitably be necessary.
The arguments to iclean will be an expression to "clean up" and a function to be Collected on.
We begin by defining an identity function symind whose eventual purpose will be to apply
symmetrization rules to an expression. We define it here so that iclean will not have to know about
the specific tensors we will be working with; it will be redefined below as we teach Mathematica about
General Relativity:
iclean[ expr_, coll_] := Block[ {ex, modex, len, res, ftab, funcs, inds, uinds},
If[ Head[ex] === Plus, len = Length[ex], (ex = {ex}; len = 1;)];
res = 0;
Do[(
Select[ modex, (Head[#1] =!= coll)&]]
/. var -> Sequence],
inds = Cases[ ftab //. Thread[ funcs -> Sequence], _Symbol];
uinds = unsun[ inds];
inds = Select[ uinds, (Count[inds, #1] > 1)&];
indi = 0;
res = res + (ex[[iind]] /. ReleaseHold[ Thread[ inds -> Hold[ind]]]);),
{iind, len}];
Now it is time to teach Mathematica something about General Relativity. We begin by defining the
variation of the square root of the determinant of the metric, using the ind function to
generate summation indices:
= det(g) ga b δ ga b
gi[a_, b_] -> gi[ Hold[ Sort[{a, b}]]] /.
G[a_, b_, c_] -> G[ Hold[ Join[ {a}, Sort[{b, c}]]]] /.
R[a_, b_] -> R[ Hold[ Sort[{a, b}]]]] /. List -> Sequence];
We can now begin to derive Einstein's Equations. We define the Einstein-Hilbert Lagrangian:
ind := Block[ {abet, num, str, index}, (
If the variable "indi" has not been used before, it is initialzed to 1. Otherwise it is incremented
and used to construct the new index symbol. The first such symbol is "aa", followed by "bb", etc. After
"zz" is generated the next index will be "aaa", and so on.
We now define a function whose purpose is to return the unsorted union of a list, removing redundant members:
abet = {"a", "b", "c", "d", "e", "f", "g", "h", "i", "j", "k", "l",
"m", "n", "o", "p", "q", "r", "s", "t", "u", "v", "w", "x", "y", "z"};
If[ Head[indi] === Symbol, indi = 1, indi = indi + 1];
unsun[ l_List] := Reap[ Sow[ 1, Flatten[l]], _, #1&][[2]];
This function is an example of "programming by unintended consequences". We eschew this programming technique
in general but borrow it here from the Mathematica Help Browser because it was too arcane to pass up
(see the Further Examples under the Reap built-in function).
symind[ expr_] := expr;
"ex" contains either the expansion of the argument expression, or if it was a simple product, a list whose
only member is the argument. The loop processes each term in the expansion.
"modex" is a modified version of "ex", where all occurences of arithmetic, derivative or
determinant functions have been eliminated; this leaves only those function invocations which we wish
to simplify. "ftab" is a table of the functions appearing in this term, with the collected function first.
"funcs" is a list of the function names, and "inds" is a list of the indices which are arguments to them.
"uinds" is used to select from inds those indices which are repeated; in this way we do not modify free
indices. Setting "indi" to zero so that the indices generated for each term by ind will start at "aa", we
rename the indices in "inds" to "aa", "bb", etc. After summing the re-indexed terms in the original expression,
we Collect on the function provided in the second argument, and further clean up indices in functions
whose arguments are symmetric (see below). Note the use of Hold and
ReleaseHold to prevent the invocation
of ind until after the Threading is complete.
(ex = Expand[ expr];
modex = ex[[iind]] /. {Times -> List, Plus -> List, Power -> List} /.
Return[ symind[ Collect[res, coll[al__]]]];)];
{der -> List, covder -> List} /. det[g_] -> 1;
ftab = Select[ unsun[ Join[
funcs = unsun[ Table[ Head[ ftab[[jind]]], {jind, Length[ftab]}]];
Select[ modex, (Head[#1] === coll)&],
(!(NumericQ[#1] || (Head[#1] === Symbol)))&];
var[ Sqrt[- det[g]]] := - Sqrt[- det[g]] g[ aind = ind, bind = ind]
This is of course based on the definition of the determinant as the exponentiation of the the trace
of the log of the metric tensor:
var[ gi[ aind, bind]] / 2;
δ det(g) = δ exp( tr( ln(g)))
We now redefine symind to sort the symmetric indices in the metric ("g"), its inverse ("gi"), the
Christoffel Symbols ("G") and the
Ricci Tensor ("R"):
= exp( tr( ln(g))) δ tr( ln(g))
symind[ expr_] := Map[ Collect[#1, Sqrt[- det[g]]]&, ReleaseHold[ expr /.
Again note the use of Hold and ReleaseHold to prevent the index sorting from
taking place until the rule substitutions have occurred.
g[a_, b_] -> g[ Hold[ Sort[{a, b}]]] /.
L := - Sqrt[- det[g]] Rs;
where "Rs" is the scalar curvature. We next find the
variation of that Lagrangian in terms of the Ricci Tensor and the inverse metric:
varL = iclean[ var[L] /. Rs -> gi[a, b] R[a, b], var] /. gi[a_, b_] R[a_, b_] -> Rs
For arbitrary variation of the (inverse) metric, the quantity in parentheses must vanish for the
Lagrangian to be stationary: this gives us the
vacuum Einstein Equations. However, the variation of
the Ricci Tensor must also vanish, and so we express it in terms of the Christoffel Symbols (by way
of the Riemann Tensor):
Sqrt[- det[g]] (1/2 Rs g[aa, bb] - R[aa, bb]) var[ gi[aa, bb]]
- Sqrt[- det[g]] gi[aa, bb] var[ R[aa, bb]]
varR = iclean[ Sqrt[- det[g]] gi[aa, bb] var[ R[aa, bb]] /.We next re-express the derivatives in terms of covariant derivatives (since the variation of the Christoffel Symbols is a tensor):R[a_, b_] -> R[a, cind = ind, b, cind] /.R[a_, b_, c_, d_] -> der[ G[d, a, c], x[b]] - der[ G[d, b, c], x[a]] +
G[find = ind, a, c] G[d, find, b] - G[find, b, c] G[d, find, a], gi]
Sqrt[- det[g]] gi[ aa, bb] (der[ var[ G[cc, aa, bb]], x[cc]] -der[ var[ G[cc, bb, cc]], x[aa]] +G[cc, cc, dd] var[ G[ dd, aa, bb]] - G[cc, bb, dd] var[ G[dd, aa, cc]] -
G[cc, aa, dd] var[ G[dd, bb, cc]] + G[cc, aa, bb] var[ G[dd, cc, dd]])
varRa = iclean[ varR /. der[ var[ G[a_, b_, c_]], x[d_]] -> covder[ var[ G[a, b, c]], x[d]] -Since the covariant derivative of the (inverse) metric is zero (we are assuming a symmetric metric), we can bring it inside the covder function:G[a, eind = ind, d] var[ G[eind, b, c]] + G[eind, c, d] var[ G[a, eind, b]] +G[eind, b, d] var[ G[a, eind, c]], gi]
(covder[ var[ G[cc, aa, bb]], x[cc]] -covder[ var[ G[cc, bb, cc]], x[aa]]) Sqrt[- det[g]] gi[aa, bb]
varRb = Expand[ varRa] /. covder[a_, b_] gi[c__] -> covder[ gi[c]a, b]These are covariant divergences of vectors which can converted to ordinary divergences by noting thatcovder[ gi[aa, bb] var[ G[cc, aa, bb]], x[cc]] Sqrt[- det[g]] -covder[ gi[aa, bb] var[ G[cc, bb, cc]], x[aa]] Sqrt[- det[g]]
(- det(g))1/2 Va; a = ((- det(g))1/2 Va), awhere the comma denotes an ordinary derivative. We examine the left hand side of this equation, expanding the Christoffel Symbol in terms of the metric and its inverse:
iclean[ Sqrt[- det[g]] (der[ vu[a], x[a]] + G[a, bind = ind, a] vu[bind]) /.and see that it is the same as the right hand side by using the variation of the determinant as an ordinary derivative:G[a_, b_, c_] -> gi[a, find = ind] * (der[ g[c, find], x[b]] +der[ g[b, find], x[c]] - der[ g[b, c], x[find]]) / 2, Sqrt]
Sqrt[- det[g]] (der[ vu[aa], x[aa]] +1/2 der[ g[aa, bb], x[cc]] gi[aa, bb] vu[cc])
iclean[ Simplify[ der[ Sqrt[- det[g]] vu[a], x[a]] /.Applying this equivalence, we see that the variation of the Ricci Tensor is an ordinary divergence, which vanishes on a manifold without boundary (or if the variations of the metric vanish on the boundary):der[ Sqrt[- det[g]], x[a]] -> - der[ det[g], x[a]] / (2 Sqrt[- det[g]]) /.der[ det[g], x[a]] -> det[ g] der[ g[aa, bb], x[a]] gi[aa, bb]], Sqrt]
Sqrt[- det[g]] (der[ vu[aa], x[aa]] +1/2 der[ g[aa, bb], x[cc]] gi[aa, bb] vu[cc])
varRb /. Sqrt[- det[g]] covder[a_, b_] -> d[ Sqrt[- det[g]]a, b]While these techniques are powerful, it is important to examine every result carefully: it is still easy to accidentally make a mess of the indices, which can result in corruption or loss of terms!d[ Sqrt[- det[g]] gi[aa, bb] var[ G[cc, aa, bb]], x[cc]] -d[ Sqrt[- det[g]] gi[aa, bb] var[ G[cc, bb, cc]], x[aa]]
It is also possible to use these methods to work with Lie Algebras. When doing so, however, you
must ensure that Mathematica does not alter the operator ordering by clearing the
Orderless attribute on the Times function:
Expand[ - z y x /. (y x) -> x y - Comm[ x, y]]
Expand[ - z y x /. (y x) -> x y - Comm[ x, y]]
The next appendix is about parallel algebraic computation.
©2005, Kenneth R. Koehler. All Rights Reserved. This document may be freely reproduced provided that this copyright notice is included.
Please send comments or suggestions to the author.
ClearAttributes[ Times, Orderless];
If this is not done, or if you accidentally reset the attribute, you will soon have operators
commuting which should not be:
- z x y + z Comm[ x, y]
SetAttributes[ Times, Orderless];
- x y z + z Comm[ x, y]