The Schwarzschild solution

In the following the Einstein field equations will be solved for the general static isotropic metric, in standard form. See Weinberg, page 175. Note: the convention for Riemann is opposite from Weinberg.

IndepVar = {t, r, θ,  φ}

The function   metric[]calculates the coefficients of the metric tensor from a given line element.   metric[] takes the line element  ds^2  and a list of coordinates as input variables. The result is the metric tensor of the underlying space.

metric[lineelement_, independentvars_List] := 
        Block[{lenindependent, differentials, diffmatrix, metricform, 
            varmetric, gh, sum, equation, rule, varhelp, zeros, zerorule}, 
          lenindependent = Length[independentvars]; 
          differentials = Dt /@ independentvars; 
          diffmatrix = Outer[Times, differentials, differentials]; 
          metricform = Array[gh, {lenindependent, lenindependent}]; 
          varmetric = Variables[metricform]; 
          If[Length[metricform] == Length[diffmatrix], sum = 0; 
            Do[Do[sum = 
                  sum + metricform[[i, 
                        j]] diffmatrix[[
                        i, j]], {j, 1, 
                  lenindependent}], {i, 1, lenindependent}], sum = 0]; 
          If[sum === 0, Return[sum], sum = sum - lineelement; 
            equation = CoefficientList[sum, differentials] == 0; 
            rule = Solve[equation, varmetric]; 
            metricform = metricform /. rule; 
            varmetric = Variables[metricform]; varhelp = {}; 
            Do[If[! (FreeQ[
                      varmetric[[i]], 
                      gh]), 
                AppendTo[varhelp, 
                  varmetric[[i]]]], {i, 
                1, Length[varmetric]}]; zeros = Table[0, {Length[varhelp]}]; 
            SubstRule[x_, y_] := x -> y; 
            zerorule = Thread[SubstRule[varhelp, zeros]]; 
            metricform = Flatten[metricform /. zerorule, 1]; 
            metricform = 
              Expand[(1/2) (metricform + Transpose[metricform])]]; 
          metricform];

In Mathematica one gets the flat cartesian metric by

\(MatrixForm[
      metric[
		\(-Dt[t]\^2\) + Dt[x]\^2 + Dt[y]\^2 + Dt[z]\^2, {t, x, y, 
          z}]]\)

The Schwarzschild metric ins standard form is

\(MatrixForm[
      metric[
		\(-B[r]\)\ Dt[t]\^2 + 
          A[r]\ Dt[r]\^2 + 
		r\^2\ Dt[θ]\^2 + 
          r\^2\ Sin[θ]\^2\ Dt[φ]\^2, {t, 
          r, θ, φ}]]\)

Equivalently one can give the line element. The line element  ds^2  in the coordinates  t ,r, th  and  phi   is

ds =  -B[r] Dt[t]^2 + A[r] Dt[r]^2 + 
        (r^2) (Dt[θ]^2 + 
              Sin[θ]^2 Dt[φ]^2)

The corresponding metric is

g = metric[ds, {t, r, θ, φ}]

ing  =  Inverse[g];

The above line element  ds^ 2  resembles the line element in Euclidean space. In the following example the radial coordinate  r  is transformed so that one can write the line element in a isotropic  form .  The corresponding line element of the metric reads

dsi =  -H[r] Dt[t]^2 + 
        J[r] (Dt[r]^2 + 
              (r^2) Dt[θ]^2 + (r^2) 
                Sin[θ]^2 Dt[φ]^2)

and the corresponding metric tensor follows by

gi = metric[dsi, {t, r, θ, φ}]

In the following we will work with the metric in "standard" form. Up to now we have only discussed the line element of the metric and its related metric tensor. To derive the field equations for the vacuum case in GR, we have to introduce other tensors. One of the essential quantities determining the field equations are the Christoffel symbols. These symbols are related to the metric tensor in a straightforward way.

Christoffel[m_, a_, b_, g_, ing_] := 
      Block[{n}, 
        Expand[(1/2) (∑_(n = 1)^(Length[g]))ing[[
                  m, n]] ((∂_(IndepVar[[b]]))g[[a, 
                        n]] + (∂_(IndepVar[[a]]))g[[b, 
                        n]] - (∂_(IndepVar[[n]]))g[[a, 
                        b]])]]

Riemann[a_, b_, c_, d_, g_, ing_] := 
      Block[{}, 
        Expand[(∂_(IndepVar[[c]]))Christoffel[a, b, d, g, 
                ing] - (∂_(IndepVar[[d]]))Christoffel[a, b, c, g, 
                ing] + (∑_(e = 1)^(Length[g]))Christoffel[e, b, d, 
                  g, ing] Christoffel[a, e, c, g, 
                  ing] - (∑_(e = 1)^(Length[g]))Christoffel[e, b, c, 
                  g, ing] Christoffel[a, e, d, g, ing]]]

Ricci[m_, q_, g_, ing_] := 
      Block[{a}, 
        Expand[(∑_(a = 1)^(Length[g]))Riemann[a, m, a, q, g, 
              ing]]]

RicciScalar[g_, ing_] := 
        Block[{}, 
          Expand[(∑_(a = 1)^(Length[g]))(∑_(b = 1)^(Length[g]))ing[[a, b]] Ricci[a, 
                    b, g, ing]]];

Einstein[m_, n_, g_, ing_] := 
      Ricci[m, n, g, ing] - 
        (1/2) RicciScalar[g, ing] g[[m, 
            n]]

Bianchi[a_, g_, ing_] := 
        Block[{}, 
          Expand[(∑_(n = 1)^(Length[g]))(∂_(IndepVar[[n]]))((∑_(m = 1)^(Length[g]))ing[[n, m]] Einstein[m, a, g, 
                          ing]) + (∑_(n = 1)^(Length[g]))(∑_(m = 1)^(Length[g]))Christoffel[n, m, n, g, 
                      ing] (∑_(l = 1)^(Length[g]))ing[[m, l]] Einstein[l, a, g, 
                          ing] - (∑_(n = 1)^(Length[g]))(∑_(m = 1)^(Length[g]))Christoffel[n, m, a, g, 
                      ing] (∑_(l = 1)^(Length[g]))ing[[m, l]] Einstein[l, n, g, 
                          ing]]];

\(Simplify[
	
      Table[Christoffel[i, j, k, g, ing], {i, 1, 4}, {j, 1, 4}, {k, 1, 
          4}]]\)

\(Simplify[
	
      Table[Ricci[i, j, g, ing], {i, 1, 4}, {j, 1, 4}]]\)

Simplify[Ricci[1, 1, g, ing]]

Simplify[Ricci[2, 2, g, ing]]

Simplify[Ricci[3, 3, g, ing]]

Simplify[Ricci[4, 4, g, ing]]

Simplify[
      Ricci[1, 1, g, ing]/B[r] + Ricci[2, 2, g, ing]/A[r] ]

Simplify[RicciScalar[g,  ing]]

Simplify[
      Table[Einstein[a, b, g, ing], {a, 1, 4}, {b, 1, 4}]]

With a cosmological constant λ the equations become

Simplify[
      Table[Einstein[a, b, g, ing]  +  λ g[[a, b]], {a, 1, 
          4}, {b, 1, 4}]]

A solution of the above (λ=0) equations is given by

DSolve[ -A[r]  + A[r]^2 + r A'[r]  == 0,  A[r],  
      r]

while for λ not zero one has

DSolve[ 
      A[r]  + (-1 + (r^2) λ) A[r]^2 - r A'[r]  == 
        0,  A[r],  r]

B[r]  =  1  -  ((2  M G)/r)

A[r]  =  (1/(1 - ((2  M G)/r)))

where M is an integration constant. One can verify that all components of Ricci vanish

Simplify[Table[Ricci[a, b, g, ing], {a, 1, 4}, {b, 1, 4}]]

Simplify[ RicciScalar[g, ing] ]

The Riemann tensor itself is singular at the origin, and on the horizon at r=2MG,

Simplify[ Riemann[1, 2, 1, 2, g, ing] ]

Simplify[ Riemann[2, 3, 2, 3, g, ing]]

Simplify[ Riemann[3, 4, 3, 4, g, ing]]

Simplify[ Riemann[1, 3, 1, 3, g, ing]]

Simplify[Riemann[1, 4, 1, 4, g, ing]]