# Functional Computation

Many of the results quoted above (especially in the section on warped product spaces), were derived using Mathematica as a rule-based environment for functional analysis. In this appendix we demonstrate these techniques by deriving Einstein's Equations from the Einstein-Hilbert Lagrangian.

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

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

We repeat the process for the variational operator (usually denoted δ):

Clear[ var];

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]}];

We further specify that the variation of a derivative is the derivative of the variation:
var[ der[ expr_, x_]] := der[ var[expr], x];

var[ covder[ expr_, x_]] := covder[ 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.
ind := Block[ {abet, num, str, index}, (
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];

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]];)];

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:
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;

iclean[ expr_, coll_] := Block[ {ex, modex, len, res, ftab, funcs, inds, uinds},

(ex = Expand[ expr];

If[ Head[ex] === Plus, len = Length[ex], (ex = {ex}; len = 1;)];

res = 0;

Do[(

modex = ex[[iind]] /. {Times -> List, Plus -> List, Power -> List} /.
{der -> List, covder -> List} /. det[g_] -> 1;
ftab = Select[ unsun[ Join[

/. var -> Sequence],

funcs = unsun[ Table[ Head[ ftab[[jind]]], {jind, Length[ftab]}]];

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}];

Return[ symind[ Collect[res, coll[al__]]]];)];
"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.

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:

var[ Sqrt[- det[g]]] := - Sqrt[- det[g]] g[ aind = ind, bind = ind]
var[ gi[ aind, bind]] / 2;
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:
δ det(g) = δ exp( tr( ln(g)))
= exp( tr( ln(g))) δ tr( ln(g))

= det(g) ga b δ ga b

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"):
symind[ expr_] := Map[ Collect[#1, Sqrt[- det[g]]]&, ReleaseHold[ expr /.
g[a_, b_] -> g[ Hold[ Sort[{a, 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];

Again note the use of Hold and ReleaseHold to prevent the index sorting from taking place until the rule substitutions have occurred.
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
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]]
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):
varR = iclean[ Sqrt[- det[g]] gi[aa, bb] var[ R[aa, bb]] /.
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]])

We next re-express the derivatives in terms of covariant derivatives (since the variation of the Christoffel Symbols is a tensor):
varRa = iclean[ varR /. der[ var[ G[a_, b_, c_]], x[d_]] -> covder[ var[ G[a, b, c]], x[d]] -
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]
Since the covariant derivative of the (inverse) metric is zero (we are assuming a symmetric metric), we can bring it inside the covder function:
varRb = Expand[ varRa] /. covder[a_, b_] gi[c__] -> covder[ gi[c]a, b]
covder[ 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]]
These are covariant divergences of vectors which can converted to ordinary divergences by noting that
(- det(g))1/2 Va; a = ((- det(g))1/2 Va), a
where 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]) /.
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])
and see that it is the same as the right hand side by using the variation of the determinant as an ordinary derivative:
iclean[ Simplify[ der[ Sqrt[- det[g]] vu[a], x[a]] /.
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])
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):
varRb /. Sqrt[- det[g]] covder[a_, b_] -> d[ Sqrt[- det[g]]a, b]
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]]
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!
ClearAttributes[ Times, Orderless];

Expand[ - z y x /. (y x) -> x y - Comm[ x, y]]

- z x y + z Comm[ x, y]
If this is not done, or if you accidentally reset the attribute, you will soon have operators commuting which should not be:
SetAttributes[ Times, Orderless];

Expand[ - z y x /. (y x) -> x y - Comm[ x, y]]

- x y z + z Comm[ x, y]

The next appendix is about parallel algebraic computation.