If you want to try it, go to the examples directory and:

ondra@syslik:~/sympy/examples$ python relativity.py

----------------------------------------

Christoffel symbols:

1/2*\nu'(r)

1/2*\nu'(r)

1/2*exp(\lambda(r))**(-1)*exp(\nu(r))*\nu'(r)

1/2*\lambda'(r)

-exp(\lambda(r))**(-1)*r

-exp(\lambda(r))**(-1)*sin(\theta)**2*r

r**(-1)

r**(-1)

-cos(\theta)*sin(\theta)

sin(\theta)**(-1)*cos(\theta)

sin(\theta)**(-1)*cos(\theta)

r**(-1)

r**(-1)

----------------------------------------

Ricci tensor:

-1/4*exp(\lambda(r))**(-1)*exp(\nu(r))*\lambda'(r)*\nu'(r)+1/2*exp(\lambda(r))**(-1)*exp(\nu(r))*(\nu'(r))'+exp(\lambda(r))**(-1)*exp(\nu(r))*r**(-1)*\nu'(r)+1/4*exp(\lambda(r))**(-1)*exp(\nu(r))*\nu'(r)**2

-1/4*\nu'(r)**2-1/2*(\nu'(r))'+r**(-1)*\lambda'(r)+1/4*\nu'(r)*\lambda'(r)

-exp(\lambda(r))**(-1)-1/2*exp(\lambda(r))**(-1)*r*\nu'(r)-(-1-sin(\theta)**(-2)*cos(\theta)**2)-sin(\theta)**(-2)*cos(\theta)**2+1/2*exp(\lambda(r))**(-1)*r*\lambda'(r)

1/2*exp(\lambda(r))**(-1)*sin(\theta)**2*r*\lambda'(r)-exp(\lambda(r))**(-1)*sin(\theta)**2+sin(\theta)**2-1/2*exp(\lambda(r))**(-1)*sin(\theta)**2*r*\nu'(r)

----------------------------------------

solve the Einstein's equations:

\lambda(r) = -log(C1+C2*r**(-1))

metric:

-(C1+C2*r**(-1)) 0 0 0

0 (C1+C2*r**(-1))**(-1) 0 0

0 0 r**2 0

0 0 0 r**2*sin(\theta)**2

It will start with a general, spherically symmetric metric:

gdd=Matrix((

(-exp(nu(r)),0,0,0),

(0, exp(lam(r)), 0, 0),

(0, 0, r**2, 0),

(0, 0, 0, r**2*sin(theta)**2)

))

notice the two unknown functions nu(r) and lam(r) and it calculates the Christoffel symbols, Riemann and Ricci tensors, then solves the differential (Einstein) equations for nu and lam and substitute it back to the metric to get the Schwarzschild metric (last matrix in the output from relativity.py).

If you look into the relativity.py, you will see, that when solving the equations, I am a little cheating in that I am assuming

nu(r) == -lam(r)

that can be shown by playing with the equations. That is so that I get just one differential equation for just one unknown function and then I plug it into the dsolve() and I am done.