```python
from sympy import *
init_session()
from IPython.display import display, Math
import util
```

    IPython console for SymPy 1.3 (Python 3.7.3-64-bit) (ground types: python)
    
    These commands were executed:
    >>> from __future__ import division
    >>> from sympy import *
    >>> x, y, z, t = symbols('x y z t')
    >>> k, m, n = symbols('k m n', integer=True)
    >>> f, g, h = symbols('f g h', cls=Function)
    >>> init_printing()
    
    Documentation can be found at http://docs.sympy.org/1.3/
    


# Christoffel Symbols
The Christoffel symbols for a given metric $g$ with components $g_\mu \nu$ in coordinates $x^\mu$ are given by
    $$\Gamma^\sigma_{\mu \nu} = 
    \frac{1}{2} g^{\sigma \rho} 
    (\frac{\partial}{\partial x^\mu} g_{\nu \rho} + 
    \frac{\partial}{\partial x^\nu} g_{\rho \mu} - \frac{\partial}{\partial x^\rho} g_{\mu \nu})$$


```python
def christoffel(g, x):
    """For a given metric/coordinates, return the Christoffel symbols as a 3-dimensional list.
    The upper index is the first index, while the lower indices are the second and third index.
    The input metric should be the metric with lower indices(upper indices are the inverse).
    """
    g_inv = g**-1
    
    n = g.shape[0]
    Gamma = [[[Integer(0) for i in range(n)] for j in range(n)] for k in range(n)]
    
    for i in range(n):
        for j in range(n):
            for k in range(j+1):
                for l in range(n):
                    Gamma[i][j][k] = Gamma[i][j][k] + Rational(1,2)*g_inv[i,l] * \
                    (g[k,l].diff(x[j]) + g[j,l].diff(x[k]) - g[j,k].diff(x[l]))
                Gamma[i][j][k] = simplify(trigsimp(expand_trig(Gamma[i][j][k])))
                Gamma[i][k][j] = Gamma[i][j][k]
    return Gamma
```

# Riemann Tensor
The Riemann tensor is given by
$$R^\rho_{\ \sigma \mu \nu} = 
\frac{\partial}{\partial x^\mu} \Gamma^\rho_{\nu \sigma} - 
\frac{\partial}{\partial x^\nu} \Gamma^\rho_{\mu \sigma} +
\Gamma^\rho_{\mu \lambda} \Gamma^\lambda_{\nu \sigma} -
\Gamma^\rho_{\nu \lambda} \Gamma^\lambda_{\mu \sigma}$$


```python
def riemann(Gamma, x):
    """Given the Christoffel symbols and coordinates, constructs the Riemann tensor as a 4-dimensional list.
    The upper index is the first index in the list, while the lower indices are the second, third and fourth
    indices in the list. """
    n = len(x)
    R = [[[[Integer(0) for i in range(n)] for j in range(n)] for k in range(n)] for l in range(n)]
    for i in range(n):
        for j in range(n):
            for k in range(n):
                for l in range(k):
                    R[i][j][k][l] = Gamma[i][l][j].diff(x[k]) - Gamma[i][k][j].diff(x[l])
                    for m in range(n):
                        R[i][j][k][l] = R[i][j][k][l] + Gamma[i][k][m]*Gamma[m][l][j] - Gamma[i][l][m]*Gamma[m][k][j]
                    R[i][j][k][l] = simplify(trigsimp(expand_trig(R[i][j][k][l])))
                    R[i][j][l][k] = - R[i][j][k][l]
    return R
```

# Ricci Tensor
The Ricci tensor is simply a contraction of the Riemann tensor:
$$R_{\mu \nu} = R^\lambda_{\ \mu \lambda \nu}$$


```python
def ricci_tensor(riem):
    """Returns the Ricci tensor as a Sympy matrix object."""
    n = len(riem)
    R = zeros(n,n)
    for i in range(n):
        for j in range(n):
            for k in range(n):
                R[i,j] = R[i,j] + riem[k][i][k][j]
            R[i,j] = simplify(trigsimp(expand_trig(R[i,j])))
    return R
```

# Ricci Scalar
The Ricci Scalar is given by the trace of the Ricci tensor:
$$R = R^\mu_{\ \mu} = g^{\mu \nu} R_{\mu \nu}$$


```python
def ricci_scalar(R, g):
    """Returns the Ricci scalar as a Sympy expression."""
    g_inv = g**-1
    return simplify(trigsimp(expand_trig((g_inv.T * R).trace())))
```

# Einstein Tensor
The Einstein tensor is given by 
$$G_{\mu \nu} = R_{\mu \nu} - \frac{1}{2} R g_{\mu \nu}$$


```python
def einstein(r_tensor, r, g):
    return simplify(trigsimp(expand_trig(r_tensor - Rational(1,2)*r*g)))
```

# Example 1


```python
x = symbols('theta phi')
theta = x[0]
a = symbols('a')
g = a**2 * diag(1, sin(theta)**2)
Gamma = christoffel(g, x)
util.display_nonzero(Gamma, x, "christoffel")
```


```python
R = riemann(Gamma, x)
```


```python
util.display_nonzero(R, x, "riemann")
```


```python
ric = ricci_tensor(R)
util.display_nonzero(ric, x, "ricci")
display(ric)
```


```python
scal = ricci_scalar(ric, g)
display(scal)
```


```python
ein = einstein(ric, scal, g)
util.display_nonzero(ein, x, "einstein")
```

    No nonzero elements of Einstein.


# Example 2


```python
x = symbols('t r theta phi')
h, w = symbols('h omega')
g = diag(-1, 1 + h*sin(w*(t-z)), 1 - h*sin(w*(t-z)), 1)
Gamma = christoffel(g, x)

# Expand to linear order in h
for i in range(4):
    for j in range(4):
        for k in range(4):
            Gamma[i][j][k] = Gamma[i][j][k].subs(h,0) + Gamma[i][j][k].diff(h).subs(h,0)*h

util.display_nonzero(Gamma, x, "christoffel")
```


```python
riem = riemann(Gamma, x)
for i in range(4):
    for j in range(4):
        for k in range(4):
            for l in range(4):
                riem[i][j][k][l] = riem[i][j][k][l].subs(h,0) + riem[i][j][k][l].diff(h).subs(h,0)*h
util.display_nonzero(riem, x, "riemann")
```


```python
ric = ricci_tensor(riem)
```


```python
util.display_nonzero(ric, x, "ricci")
```


```python
scal = ricci_scalar(ric, g)
```


```python
scal
```


```python
x = symbols("r theta")
r = x[0]
g = diag(1, r**2)
Gamma = christoffel(g, x)
riem = riemann(Gamma, x)
ric = ricci_tensor(riem)
scal = ricci_scalar(ric, g)
util.display_nonzero(Gamma, x, "christoffel")
util.display_nonzero(riem, x, "riemann")
util.display_nonzero(ric, x, "ricci")
display(Math(r'R = '), scal)
```

bf70839dd82e1fb30a2c6abee0085a2d5b9909a8
ipynb
Jupyter Notebook
gr_tensors.ipynb
rahulmalayappan/general-relativity-tensors
bae2ccde50875d34e40f8a3ec33cb1bc0324c943
MIT
gr_tensors.ipynb
rahulmalayappan/general-relativity-tensors
bae2ccde50875d34e40f8a3ec33cb1bc0324c943
MIT
gr_tensors.ipynb
rahulmalayappan/general-relativity-tensors
bae2ccde50875d34e40f8a3ec33cb1bc0324c943
MIT
Qwen/Qwen-72B
1. YES
2. YES