Curvature and the geodesic equations

The code in this notebook is based on James Hartle’s Mathematica programs written as companions to his book “Gravity: An Introduction to Einstein’s General Relativity”, available from the Book’s site:
http://web.physics.ucsb.edu/~gravitybook/mathematica.html
The original program was written by Leonard Parker, University of Wisconsin, Milwaukee (see the Acknowledgement section in the end of this notebook). The code in this notebook is a modification of the original program, adapted for the needs of the course “General Relativity and Cosmology” offered at the National Technical University of Athens, by Konstantinos Anagnostopoulos ( http://physics.ntua.gr/konstant ).The site of the course for the Spring 2023 can be found at  http://physics.ntua.gr/konstant/GR .

This is the Mathematica notebook Curvature and the Einstein Equation available from the book website.  From a given metric (Subscript[g, αβ]) , it computes the components of the following: the inverse metric, (g^(λσ)), the Christoffel symbols or affine connection,

(Subscript[(Γ^λ), μν])=(1/2)(g^(λσ))((∂_μ) (Subscript[g, σν])+(∂_ν) (Subscript[g, σμ])-(∂_σ) (Subscript[g, μν])),

( (∂_(α  ))stands for the partial derivative ∂/∂(x^α)), the Riemann tensor,

(Subscript[(R^λ), μνσ])=(∂_ν) (Subscript[(Γ^λ), μσ])-(∂_σ) (Subscript[(Γ^λ), μν])+(Subscript[(Γ^λ), ην])(Subscript[(Γ^η), μσ])-(Subscript[(Γ^λ), ησ])(Subscript[(Γ^η), μν]),

the Ricci tensor

(Subscript[R, μν])=(Subscript[(R^λ), μλν]),

the scalar curvature,

R=(g^(μν))(Subscript[R, μν]),

and the Einstein tensor,

(Subscript[G, μν])=(Subscript[R, μν])-(1/2)(Subscript[g, μν])R.

We also compute:
(Subscript[R, μνρσ]), (R^(μνρσ)), (R^2)=(Subscript[R, μνρσ]) (R^(μνρσ)), the Weyl Tensor
 (Subscript[ℭ, ρσμν])=(Subscript[R, ρσμν])-(2/(n-2)){(Subscript[g, ρ[μ])(Subscript[R, ν]σ])-(Subscript[g, σ[μ])(Subscript[R, ν]ρ])}+(2/((n-1)(n-2)))(Subscript[g, ρ[μ])(Subscript[g, ν]σ])  R, 
 and the geodesic equations.

You must input the covariant components of the metric tensor (Subscript[g, μν]) by editing the relevant input line in this Mathematica notebook. You may also wish to change the names of the coordinates. Only the nonzero components of the above quantities are displayed as the output. All the components computed are in the coordinate basis in which the metric was specified.

Initialization:

First clear any values that may already have been assigned to the names of the various objects to be calculated. The names of the coordinates that you will use are also cleared. The  name of the coordinates must be defined (the list coord), and then the dimension of the manifold n is set.

Note:
It is important not to use the symbols, i, j, k, l, s, R, G, ℭ, τ, Γ, or n as constants or coordinates in the metric that you specify above. The reason is that the first five of  those symbols are used as summation or table indices in the calculations done below, and n is the dimension of the space. The rest are also used for displaying results in the program.

Clear[coord, metric,inversemetric, affine, riemann, lriemann, uriemann, ricci, scalar, einstein,weyl, geodesic,R,G,τ, i,j,k,l,s];
Clear[r,θ,ϕ,t,χ,a,m]; 

(*----------------------------------------------------------------------------*)
(* This is what you need to set + the metric below: *)
coord = {t,χ,θ,ϕ};
n          =Length[coord];
(*----------------------------------------------------------------------------*)
Print["The Manifold has dimension   n= ",n,"\nCoordinate system:          ",coord]

Defining the metric:   The Friedmann metric

Input the metric as a list of lists, i.e., as a matrix. You can input the components of any metric here, but you must specify them as explicit functions of the coordinates.

Compute for the metric:
(ds^2)=-(dt^2)+((a(t))^2)[ ((dχ)^2)+(sin^2)(χ) (((dθ)^2)+(sin^2)(θ)((dϕ)^2)) ]

(*----------------------------------------------------------------------------*)
metric={
{-1,0      ,0                         ,                                        0},
{0,((a[t])^2),0                         ,                                        0},
{0,0        ,((a[t])^2) ((Sin[χ])^2),                                         0},
{0,0       ,0                          ,((a[t])^2) ((Sin[χ])^2) ((Sin[θ])^2)}};
(*----------------------------------------------------------------------------*)
metric//MatrixForm

Calculating the inverse metric:

The inverse metric is obtained through matrix inversion.

inversemetric=FullSimplify[Inverse[metric]];
inversemetric//MatrixForm

Let’s see how we can compute the derivatives:

metric[[2,2]]

coord[[1]]

D[metric[[2,2]],t]

Use the list of coordinates to specify the variable w.r.t. which we differentiate:

D[metric[[2,2]],coord[[1]]]

Table[D[metric[[i,i]],coord[[1]]],{i,n}]

Table[D[metric[[i,i]],coord[[j]]],{i,n},{j,n}]

(∂_t)metric[[2,2]]

(∂_(coord[[1]]))metric[[2,2]]

Calculating the Christoffel symbols:

The calculation of the components of the Christoffel symbols is done by transcribing the definition given earlier into the notation of Mathematica and using the Mathematica functions D for taking partial derivatives, Sum for summing over repeated indices, Table for forming a list of components, and Simplify for simplifying the result.

The delayed definition performs a calculation on demand and caches the result:
E.g.: pi2:=pi2=N[Pi^2,50];
fib[1]=fib[2]=1;
fib[n_]:=fib[n]=fib[n-1]+fib[n-2]
fib[6] -> 8
Definition[fib]
({{fib[1]=1}, {fib[2]=1}, {fib[3]=2}, {fib[4]=3}, {fib[5]=5
fib[6]=8}, {fib[n_]:=fib[n]=fib[n-1]+fib[n-2]}})

affine:=affine=FullSimplify[Table[
(1/2)*Sum[
(*              (g^is)    (  (∂_k)(Subscript[g, sj])                   +   (∂_j)(Subscript[g, sk])                  -     (∂_s)(Subscript[g, jk]))             *)
(inversemetric[[i,s]])*(D[metric[[s,j]],coord[[k]] ]+D[metric[[s,k]],coord[[j]] ]-D[metric[[j,k]],coord[[s]] ]),
{s,1,n}],
{i,1,n},{j,1,n},{k,1,n}] ];

Now let’s see how to display the nontrivial results. First let’s look what affine consists of:

affine

We can look at (Subscript[(Γ^μ), νρ]) as matrices (Subscript[((Γ^μ)), νρ]):

Table[affine[[i]]//MatrixForm,{i,1,n}]

Many zeroes... We also need to associate the nonzero values with the indices μ, ν, ρ. Let’s make pairs of values and (Subscript[(Γ^μ), νρ]) symbols:
We use the functions Superscript[Γ,i]->(Γ^i), and Subscript[A,j,k]->(Subscript[A, j,k]). The function UnsameQ[<expression>,0] evaluates to True, when <expression> is not the same as 0. If[<condition>,<value>] returns <value> when <condition> is True and Null when it is False:

listaffine:=Table[
If[
UnsameQ[affine[[i,j,k]],0],
{Subscript[Superscript[Γ,i],j,k],affine[[i,j,k]]}
  ] ,
{i,1,n},{j,1,n},{k,1,j}];

Now let’s see what listaffine consists of:

listaffine

Now we want to get rid of the Nulls: First Flatten the list:

Flatten[listaffine]

Then get rid of the Nulls, using the DelteCases function:

DeleteCases[Flatten[listaffine],Null]

Then, reconstruct pairs of {Symbol,Values} using the Partition function:

Partition[DeleteCases[Flatten[listaffine],Null],2]

Yeah, we are ready, now let’s make a pretty display:

Partition[DeleteCases[Flatten[listaffine],Null],2]//TableForm

Or use more spacing for more comfortable reading: TableSpacing -> {yspaces, xspaces} adds spaces in rows+columns

TableForm[Partition[DeleteCases[Flatten[listaffine],Null],2],TableSpacing->{4,8}]

Calculating and displaying the Riemann tensor:

The components of the Riemann tensor, (Subscript[(R^λ), μνσ]), are calculated using the definition given above.
(Subscript[(R^i), jkl])=(∂_k)(Subscript[(Γ^i), lj])-(∂_l)(Subscript[(Γ^i), kj])+(Subscript[(Γ^i), km]) (Subscript[(Γ^m), lj]) - (Subscript[(Γ^i), lm]) (Subscript[(Γ^m), kj])=(Subscript[(Γ^i), lj ; k])-(Subscript[(Γ^i), kj ; l])+(Subscript[(Γ^i), ks]) (Subscript[(Γ^s), lj]) - (Subscript[(Γ^i), ls]) (Subscript[(Γ^s), kj])

riemann:=riemann=FullSimplify[Table[
(* (Subscript[(R^i), jkl])=    (∂_k)(Subscript[(Γ^i), lj])                   -         (∂_l)(Subscript[(Γ^i), kj])*)
D[       affine[[i,l,j]],  coord[[k]]           ]-D[affine[[i,k,j]], coord[[l]]          ]+
(*         (Subscript[(Γ^i), ks])             (Subscript[(Γ^s), lj])     -         (Subscript[(Γ^i), ls])            (Subscript[(Γ^s), kj])      *)
Sum[affine[[i,k,s]]affine[[s,l,j]] -    affine[[i,l,s]] affine[[s,k,j]] ,
{s,1,n}],
{i,1,n},{j,1,n},{k,1,n},{l,1,n}] ];
listriemann:=Table[
If[
UnsameQ[riemann[[i,j,k,l]],0],
{Subscript[Superscript[R,i],j,k,l],riemann[[i,j,k,l]]}
] ,
{i,1,n},{j,1,n},{k,1,n},{l,1,k-1}];
TableForm[Partition[DeleteCases[Flatten[listriemann],Null],2],TableSpacing->{2,2}]

Table[riemann[[i,j]]//MatrixForm,{i,n},{j,n}]

Check that: (Subscript[(R^μ), ν ρσ])=-(Subscript[(R^μ), ν σρ])

Table[riemann[[i,j]]//MatrixForm,{i,n},{j,n}]//MatrixForm

riemann[[i,j]] gives the matrix (Subscript[((Subscript[(R^i), j])), kl])
e.g. {{a,b,c},{d,e,f},{g,h,i}}[[2,3]] -> f  (2->row 3->column)

riemann[[1,2,2,1]]

Now we calculate the (Subscript[R, μνρλ])=(Subscript[g, μσ])(Subscript[(R^σ), νρλ])

lriemann:=lriemann=FullSimplify[Table[
Sum[metric[[i,ii]] riemann[[ii,j,k,l]],{ii,1,n}],
{i,1,n},{j,1,n},{k,1,n},{l,1,n}]];
listlriemann:=Table[
If[
UnsameQ[lriemann[[i,j,k,l]],0],
{Subscript[R,i,j,k,l],lriemann[[i,j,k,l]]}
] ,{i,1,n},{j,1,i-1},{k,1,n},{l,1,k-1}];
TableForm[Partition[DeleteCases[Flatten[listlriemann],Null],2],TableSpacing->{2,2}]

Check the following antisymmetry properties: (Subscript[R, μν ρσ])=-(Subscript[R, μν σρ])=-(Subscript[R, ⋁μ ρσ])

Table[lriemann[[i,j]]//MatrixForm,{i,n},{j,n}]//MatrixForm

Now calculate contravariant Riemann: (R^(μνρλ))=(g^(⋁α))(g^(ρβ))(g^(λγ))(Subscript[(R^μ), αβγ])

uriemann:=uriemann=FullSimplify[Table[
Sum[
inversemetric[[j,jj]] inversemetric[[k,kk]] inversemetric[[l,ll]] riemann[[i,jj,kk,ll]],{jj,1,n},{kk,1,n},{ll,1,n}
],
{i,1,n},{j,1,n},{k,1,n},{l,1,n}]];
listuriemann:=Table[
If[
UnsameQ[uriemann[[i,j,k,l]],0],{Superscript[Superscript[Superscript[Superscript[R,i],j],k],l],uriemann[[i,j,k,l]]}
] ,{i,1,n},{j,1,i-1},{k,1,n},{l,1,k-1}];
TableForm[Partition[DeleteCases[Flatten[listuriemann],Null],2],TableSpacing->{2,2}]

Calculating and displaying the (R^2) scalar:

Now calculate (R^2)=(Subscript[R, μνλρ])(R^(μνλρ))

r2=FullSimplify[Sum[lriemann[[i,j,k,l]]uriemann[[i,j,k,l]],{i,1,n},{j,1,n},{k,1,n},{l,1,n}]]

Calculating and displaying the Ricci tensor:

The Ricci tensor (Subscript[R, μν]) was defined by summing the first and third indices of the Riemann tensor (which has the first index already raised).

ricci:=ricci=FullSimplify[Table[
Sum[
riemann[[i,j,i,l]],
{i,1,n}
],{j,1,n},{l,1,n}] ];
listricci:=Table[
If[
UnsameQ[ricci[[j,l]],0],
{Subscript[R,j,l],ricci[[j,l]]}
] ,{j,1,n},{l,1,j}];
TableForm[Partition[DeleteCases[Flatten[listricci],Null],2],TableSpacing->{2,2}]

ricci//MatrixForm

A vanishing table (as with the Schwarzschild metric example)  means that the vacuum Einstein equation is satisfied.

Calculating the scalar curvature:

The scalar curvature R is calculated using the inverse metric and the Ricci tensor. The result is displayed in the output line.

scalar=FullSimplify[Sum[inversemetric[[i,j]]ricci[[i,j]],{i,1,n},{j,1,n}] ]

Calculating the Einstein tensor:

The Einstein tensor,  (Subscript[G, μν])=(Subscript[R, μν])-(1/2)(Subscript[g, μν])R, is found from the tensors already calculated.

einstein:=einstein=FullSimplify[ricci-(1/2)scalar*metric];
listeinstein:=Table[
If[
UnsameQ[einstein[[j,l]],0],
{Subscript[G,j,l],einstein[[j,l]]}
] ,{j,1,n},{l,1,j}];
TableForm[Partition[DeleteCases[Flatten[listeinstein],Null],2],TableSpacing->{2,2}]

einstein //MatrixForm

A vanishing table means that the vacuum Einstein equation is satisfied!

Calculating the Weyl tensor:

The Weyl tensor: 
(Subscript[ℭ, ρσμν])=(Subscript[R, ρσμν])-(2/(n-2)){(Subscript[g, ρ[μ])(Subscript[R, ν]σ])-(Subscript[g, σ[μ])(Subscript[R, ν]ρ])}                     +(2/((n-1)(n-2)))(Subscript[g, ρ[μ])(Subscript[g, ν]σ])                R
(Subscript[ℭ, ijkl])    =(Subscript[R, ijkl])    -(1/(n-2)){(Subscript[g, ik])(Subscript[R, lj])-(Subscript[g, il])(Subscript[R, kj])-(Subscript[g, jk])(Subscript[R, li])+(Subscript[g, jl])(Subscript[R, ki])}+(1/((n-1)(n-2))){(Subscript[g, ik])(Subscript[g, lj])-(Subscript[g, il])(Subscript[g, kj])}R

weyl:=weyl=FullSimplify[Table[
If[n>3,
lriemann[[i,j,k,l]]
-         (1/(n-2))           (metric[[i,k]] ricci[[l,j]] -metric[[i,l]] ricci[[k,j]]-metric[[j,k]] ricci[[l,i]]+metric[[j,l]] ricci[[k,i]])
+(1/((n-1)(n-2)))(metric[[i,k]] metric[[l,j]]-metric[[i,l]] metric[[k,j]]) scalar
(*else, if n<= 3 return 0:*),0],
{i,1,n},{j,1,n},{k,1,n},{l,1,n}]];
listweyl:=Table[
If[
UnsameQ[weyl[[i,j,k,l]],0],
{Subscript[ℭ,i,j,k,l],weyl[[i,j,k,l]]}
] ,{i,1,n},{j,1,i-1},{k,1,n},{l,1,k-1}];
TableForm[Partition[DeleteCases[Flatten[listweyl],Null],2],TableSpacing->{2,2}]

Calculating the geodesic equations:

The geodesic equations are calculated by asking  Mathematica to carry out the sum -(Subscript[(Γ^α), βγ])(u^β)(u^γ),  where (u^α) are the components of the four-velocity.  (This gives the derivative of (u^α) with respect to proper time  τ.  (This is replaced by  s  if the geodesics are spacelike.)

geodesic:=geodesic=Simplify[Table[
-Sum[affine[[i,j,k]]u[j]u[k],{j,1,n},
{k,1,n}],{i,1,n}]];
listgeodesic:=Table[{"d/dτ"ToString[u[i]],"=",geodesic[[i]]},{i,1,n}];
TableForm[listgeodesic,TableSpacing->{2}]

And here they are displayed as differential equations for the coordinates. τ is the affine parameter.

geodesic:=geodesic=Simplify[Table[-Sum[affine[[i,j,k]]u[j]u[k],{j,1,n},{k,1,n}],{i,1,n}]];
subst = Table[u[i]-> Subscript[coord[[i]],τ],{i,1,n}];
nlistgeodesic:=Table[{Subscript[coord[[i]],ττ],"+",-geodesic[[i]]/.subst,"= 0"},{i,1,n}];
TableForm[nlistgeodesic,TableSpacing->{2}]

The Sphere (S^2)

We use coordinates (θ,ϕ), and the metric:

(ds^2)=((dθ)^2)+(sin^2)(θ) ((dϕ)^2)

Clear[coord, metric,inversemetric, affine, riemann, lriemann, uriemann, ricci, scalar, einstein,weyl, geodesic,R,G,τ, i,j,k,l,s,m];
Clear[r,θ,ϕ,t,χ,a]; 

(*----------------------------------------------------------------------------*)
(* This is what you need to set: *)
coord = {θ,ϕ};
n          =Length[coord];
metric={
{1,0             },
{0,((Sin[θ])^2)}
};
(*----------------------------------------------------------------------------*)

inversemetric=FullSimplify[Inverse[metric]];
Print["-----------------------------------------------------------------------"];
Print["The Manifold has dimension   n= ",n,"\nCoordinate system:          ",coord];
Print["-----------------------------------------------------------------------"];
Print["(Subscript[g, μν])=",metric//MatrixForm];
Print["(g^(μν))=",inversemetric//MatrixForm];
affine:=affine=FullSimplify[Table[
(1/2)*Sum[
(*              (g^is) ((∂_k)(Subscript[g, sj])+(∂_j)(Subscript[g, sk])-(∂_s)(Subscript[g, jk]))           *)
(inversemetric[[i,s]])*
(D[metric[[s,j]],coord[[k]] ]+D[metric[[s,k]],coord[[j]] ]-D[metric[[j,k]],coord[[s]] ]),
{s,1,n}],
{i,1,n},{j,1,n},{k,1,n}] ];
(*The non zero Christoffel symbols are computed and selected below: *)
listaffine:=Table[
If[
UnsameQ[affine[[i,j,k]],0],
{Subscript[Superscript[Γ,i],j,k],affine[[i,j,k]]}
  ] ,
{i,1,n},{j,1,n},{k,1,j}];
Print["-----------------------------------------------------------------------"];
Print["Christoffel Symbols:"];
Print[TableForm[Partition[DeleteCases[Flatten[listaffine],Null],2],TableSpacing->{2,2}]];
riemann:=riemann=FullSimplify[Table[
(* (Subscript[(R^i), jkl])=    (∂_k)(Subscript[(Γ^i), lj])                   -         (∂_l)(Subscript[(Γ^i), kj])*)
D[       affine[[i,l,j]],  coord[[k]]           ]-D[affine[[i,k,j]], coord[[l]]          ]+
(*         (Subscript[(Γ^i), ks])             (Subscript[(Γ^s), lj])     -         (Subscript[(Γ^i), ls])            (Subscript[(Γ^s), kj])      *)
Sum[affine[[i,k,s]]affine[[s,l,j]] -    affine[[i,l,s]] affine[[s,k,j]] ,
{s,1,n}],
{i,1,n},{j,1,n},{k,1,n},{l,1,n}] ];
listriemann:=Table[
If[
UnsameQ[riemann[[i,j,k,l]],0],
{Subscript[Superscript[R,i],j,k,l],riemann[[i,j,k,l]]}
] ,
{i,1,n},{j,1,n},{k,1,n},{l,1,k-1}];
Print["-----------------------------------------------------------------------"];
Print["Riemann Tensor:"];
Print[TableForm[Partition[DeleteCases[Flatten[listriemann],Null],2],TableSpacing->{2,2}]];
lriemann:=lriemann=FullSimplify[Table[
Sum[metric[[i,ii]] riemann[[ii,j,k,l]],{ii,1,n}],
{i,1,n},{j,1,n},{k,1,n},{l,1,n}]];
listlriemann:=Table[
If[
UnsameQ[lriemann[[i,j,k,l]],0],
{Subscript[R,i,j,k,l],lriemann[[i,j,k,l]]}
] ,{i,1,n},{j,1,i-1},{k,1,n},{l,1,k-1}];
Print["-----------------------------------------------------------------------"];
Print["Contravariant Riemann Tensor:"];
Print[TableForm[Partition[DeleteCases[Flatten[listlriemann],Null],2],TableSpacing->{2,2}]];
uriemann:=uriemann=FullSimplify[Table[
Sum[
inversemetric[[j,jj]] inversemetric[[k,kk]] inversemetric[[l,ll]] riemann[[i,jj,kk,ll]],{jj,1,n},{kk,1,n},{ll,1,n}
],
{i,1,n},{j,1,n},{k,1,n},{l,1,n}]];
listuriemann:=Table[
If[
UnsameQ[uriemann[[i,j,k,l]],0],{Superscript[Superscript[Superscript[Superscript[R,i],j],k],l],uriemann[[i,j,k,l]]}
] ,{i,1,n},{j,1,i-1},{k,1,n},{l,1,k-1}];
Print["-----------------------------------------------------------------------"];
Print["Covariant Riemann Tensor:"];
Print[TableForm[Partition[DeleteCases[Flatten[listuriemann],Null],2],TableSpacing->{2,2}]];
r2=FullSimplify[Sum[lriemann[[i,j,k,l]]uriemann[[i,j,k,l]],{i,1,n},{j,1,n},{k,1,n},{l,1,n}]];
Print["-----------------------------------------------------------------------"];
Print["(R^2)= ",r2];
ricci:=ricci=FullSimplify[Table[
Sum[
riemann[[i,j,i,l]],
{i,1,n}
],{j,1,n},{l,1,n}] ];
listricci:=Table[
If[
UnsameQ[ricci[[j,l]],0],
{Subscript[R,j,l],ricci[[j,l]]}
] ,{j,1,n},{l,1,j}];
Print["-----------------------------------------------------------------------"];
Print["Ricci Tensor:"];
Print[TableForm[Partition[DeleteCases[Flatten[listricci],Null],2],TableSpacing->{2,2}]];
scalar=FullSimplify[Sum[inversemetric[[i,j]]ricci[[i,j]],{i,1,n},{j,1,n}] ];
Print["-----------------------------------------------------------------------"];
Print["Curvature Scalar:"];
Print["R= ",scalar];
einstein:=einstein=FullSimplify[ricci-(1/2)scalar*metric];
listeinstein:=Table[
If[
UnsameQ[einstein[[j,l]],0],
{Subscript[G,j,l],einstein[[j,l]]}
] ,{j,1,n},{l,1,j}];
Print["-----------------------------------------------------------------------"];
Print["Einstein Tensor:"];
Print[TableForm[Partition[DeleteCases[Flatten[listeinstein],Null],2],TableSpacing->{2,2}]];
weyl:=weyl=FullSimplify[Table[
If[n>3,
lriemann[[i,j,k,l]]
-         (1/(n-2))           (metric[[i,k]] ricci[[l,j]] -metric[[i,l]] ricci[[k,j]]-metric[[j,k]] ricci[[l,i]]+metric[[j,l]] ricci[[k,i]])
+(1/((n-1)(n-2)))(metric[[i,k]] metric[[l,j]]-metric[[i,l]] metric[[k,j]]) scalar
(*else, if n<= 3 return 0:*),0],
{i,1,n},{j,1,n},{k,1,n},{l,1,n}]];
listweyl:=Table[
If[
UnsameQ[weyl[[i,j,k,l]],0],
{Subscript[ℭ,i,j,k,l],weyl[[i,j,k,l]]}
] ,{i,1,n},{j,1,i-1},{k,1,n},{l,1,k-1}];
Print["-----------------------------------------------------------------------"];
Print["Weyl Tensor:"];
Print[TableForm[Partition[DeleteCases[Flatten[listweyl],Null],2],TableSpacing->{2,2}]];
geodesic:=geodesic=Simplify[Table[-Sum[affine[[i,j,k]]u[j]u[k],{j,1,n},{k,1,n}],{i,1,n}]];
subst = Table[u[i]-> Subscript[coord[[i]],τ],{i,1,n}];
nlistgeodesic:=Table[{Subscript[coord[[i]],ττ],"+",-geodesic[[i]]/.subst,"= 0"},{i,1,n}];
Print["-----------------------------------------------------------------------"];
Print["Geodesic Equations:"];
Print[TableForm[nlistgeodesic,TableSpacing->{2}]];

The Schwarzschild Metric

We consider the Schwarzschild metric outside the horizon. We use (t,r,θ,ϕ) coordinates, and the line element:

(ds^2)=-(1-((2m)/r)) (dt^2)+(((1-((2m)/r)))^(-1)) (dr^2) + (r^2)(((dθ)^2)+(sin^2)(θ) ((dϕ)^2))

Clear[coord, metric,inversemetric, affine, riemann, lriemann, uriemann, ricci, scalar, einstein,weyl, geodesic,R,G,τ, i,j,k,l,s,m];
Clear[r,θ,ϕ,t,χ,a]; 

(*----------------------------------------------------------------------------*)
(* This is what you need to set: *)
coord = {t,r,θ,ϕ};
n          =Length[coord];
metric={
{-(1- ((2m)/r)),0           ,0       ,                 0},
{0                   ,(1/(1-((2m)/r))),0      ,                   0},
{0                   ,0           ,(r^(2     )),                    0},
{0                   ,0           ,0       ,(r^2) ((Sin[θ])^2)}};
(*----------------------------------------------------------------------------*)

inversemetric=FullSimplify[Inverse[metric]];
Print["-----------------------------------------------------------------------"];
Print["The Manifold has dimension   n= ",n,"\nCoordinate system:          ",coord];
Print["-----------------------------------------------------------------------"];
Print["(Subscript[g, μν])=",metric//MatrixForm];
Print["(g^(μν))=",inversemetric//MatrixForm];
affine:=affine=FullSimplify[Table[
(1/2)*Sum[
(*              (g^is) ((∂_k)(Subscript[g, sj])+(∂_j)(Subscript[g, sk])-(∂_s)(Subscript[g, jk]))           *)
(inversemetric[[i,s]])*
(D[metric[[s,j]],coord[[k]] ]+D[metric[[s,k]],coord[[j]] ]-D[metric[[j,k]],coord[[s]] ]),
{s,1,n}],
{i,1,n},{j,1,n},{k,1,n}] ];
(*The non zero Christoffel symbols are computed and selected below: *)
listaffine:=Table[
If[
UnsameQ[affine[[i,j,k]],0],
{Subscript[Superscript[Γ,i],j,k],affine[[i,j,k]]}
  ] ,
{i,1,n},{j,1,n},{k,1,j}];
Print["-----------------------------------------------------------------------"];
Print["Christoffel Symbols:"];
Print[TableForm[Partition[DeleteCases[Flatten[listaffine],Null],2],TableSpacing->{2,2}]];
riemann:=riemann=FullSimplify[Table[
(* (Subscript[(R^i), jkl])=    (∂_k)(Subscript[(Γ^i), lj])                   -         (∂_l)(Subscript[(Γ^i), kj])*)
D[       affine[[i,l,j]],  coord[[k]]           ]-D[affine[[i,k,j]], coord[[l]]          ]+
(*         (Subscript[(Γ^i), ks])             (Subscript[(Γ^s), lj])     -         (Subscript[(Γ^i), ls])            (Subscript[(Γ^s), kj])      *)
Sum[affine[[i,k,s]]affine[[s,l,j]] -    affine[[i,l,s]] affine[[s,k,j]] ,
{s,1,n}],
{i,1,n},{j,1,n},{k,1,n},{l,1,n}] ];
listriemann:=Table[
If[
UnsameQ[riemann[[i,j,k,l]],0],
{Subscript[Superscript[R,i],j,k,l],riemann[[i,j,k,l]]}
] ,
{i,1,n},{j,1,n},{k,1,n},{l,1,k-1}];
Print["-----------------------------------------------------------------------"];
Print["Riemann Tensor:"];
Print[TableForm[Partition[DeleteCases[Flatten[listriemann],Null],2],TableSpacing->{2,2}]];
lriemann:=lriemann=FullSimplify[Table[
Sum[metric[[i,ii]] riemann[[ii,j,k,l]],{ii,1,n}],
{i,1,n},{j,1,n},{k,1,n},{l,1,n}]];
listlriemann:=Table[
If[
UnsameQ[lriemann[[i,j,k,l]],0],
{Subscript[R,i,j,k,l],lriemann[[i,j,k,l]]}
] ,{i,1,n},{j,1,i-1},{k,1,n},{l,1,k-1}];
Print["-----------------------------------------------------------------------"];
Print["Contravariant Riemann Tensor:"];
Print[TableForm[Partition[DeleteCases[Flatten[listlriemann],Null],2],TableSpacing->{2,2}]];
uriemann:=uriemann=FullSimplify[Table[
Sum[
inversemetric[[j,jj]] inversemetric[[k,kk]] inversemetric[[l,ll]] riemann[[i,jj,kk,ll]],{jj,1,n},{kk,1,n},{ll,1,n}
],
{i,1,n},{j,1,n},{k,1,n},{l,1,n}]];
listuriemann:=Table[
If[
UnsameQ[uriemann[[i,j,k,l]],0],{Superscript[Superscript[Superscript[Superscript[R,i],j],k],l],uriemann[[i,j,k,l]]}
] ,{i,1,n},{j,1,i-1},{k,1,n},{l,1,k-1}];
Print["-----------------------------------------------------------------------"];
Print["Covariant Riemann Tensor:"];
Print[TableForm[Partition[DeleteCases[Flatten[listuriemann],Null],2],TableSpacing->{2,2}]];
r2=FullSimplify[Sum[lriemann[[i,j,k,l]]uriemann[[i,j,k,l]],{i,1,n},{j,1,n},{k,1,n},{l,1,n}]];
Print["-----------------------------------------------------------------------"];
Print["(R^2)= ",r2];
ricci:=ricci=FullSimplify[Table[
Sum[
riemann[[i,j,i,l]],
{i,1,n}
],{j,1,n},{l,1,n}] ];
listricci:=Table[
If[
UnsameQ[ricci[[j,l]],0],
{Subscript[R,j,l],ricci[[j,l]]}
] ,{j,1,n},{l,1,j}];
Print["-----------------------------------------------------------------------"];
Print["Ricci Tensor:"];
Print[TableForm[Partition[DeleteCases[Flatten[listricci],Null],2],TableSpacing->{2,2}]];
scalar=FullSimplify[Sum[inversemetric[[i,j]]ricci[[i,j]],{i,1,n},{j,1,n}] ];
Print["-----------------------------------------------------------------------"];
Print["Curvature Scalar:"];
Print["R= ",scalar];
einstein:=einstein=FullSimplify[ricci-(1/2)scalar*metric];
listeinstein:=Table[
If[
UnsameQ[einstein[[j,l]],0],
{Subscript[G,j,l],einstein[[j,l]]}
] ,{j,1,n},{l,1,j}];
Print["-----------------------------------------------------------------------"];
Print["Einstein Tensor:"];
Print[TableForm[Partition[DeleteCases[Flatten[listeinstein],Null],2],TableSpacing->{2,2}]];
weyl:=weyl=FullSimplify[Table[
If[n>3,
lriemann[[i,j,k,l]]
-         (1/(n-2))           (metric[[i,k]] ricci[[l,j]] -metric[[i,l]] ricci[[k,j]]-metric[[j,k]] ricci[[l,i]]+metric[[j,l]] ricci[[k,i]])
+(1/((n-1)(n-2)))(metric[[i,k]] metric[[l,j]]-metric[[i,l]] metric[[k,j]]) scalar
(*else, if n<= 3 return 0:*),0],
{i,1,n},{j,1,n},{k,1,n},{l,1,n}]];
listweyl:=Table[
If[
UnsameQ[weyl[[i,j,k,l]],0],
{Subscript[ℭ,i,j,k,l],weyl[[i,j,k,l]]}
] ,{i,1,n},{j,1,i-1},{k,1,n},{l,1,k-1}];
Print["-----------------------------------------------------------------------"];
Print["Weyl Tensor:"];
Print[TableForm[Partition[DeleteCases[Flatten[listweyl],Null],2],TableSpacing->{2,2}]];
geodesic:=geodesic=Simplify[Table[-Sum[affine[[i,j,k]]u[j]u[k],{j,1,n},{k,1,n}],{i,1,n}]];
subst = Table[u[i]-> Subscript[coord[[i]],τ],{i,1,n}];
nlistgeodesic:=Table[{Subscript[coord[[i]],ττ],"+",-geodesic[[i]]/.subst,"= 0"},{i,1,n}];
Print["-----------------------------------------------------------------------"];
Print["Geodesic Equations:"];
Print[TableForm[nlistgeodesic,TableSpacing->{2}]];

We can see the following symmetries of the Weyl tensor in this display: (Subscript[C, μν ρσ])=-(Subscript[C, μν σρ])=-(Subscript[C, νμ ρσ])

Table[weyl[[i,j]]//MatrixForm,{i,n},{j,n}]//MatrixForm

Acknowledgment

The original program was kindly written by Leonard Parker, University of Wisconsin, Milwaukee especially  for this text. 

Several additions by Konstantinos Anagnostopoulos, National Tech U. Athens, http://physics.ntua.gr/konstant, konstant@mail.ntua.gr, 2023