I have covered several different software packages for doing scientific computation in Linux Journal, but I haven't spent as much time describing available libraries and the kind of work that can be done with those libraries. So in this article, I take a look at SymPy, the Python module that allows you to do symbolic mathematics, and specifically at a utility named GraviPy that is built upon SymPy. With this utility, you can do all sorts of computations that you'll want to do if you're researching general relativity.

You can install both modules using pip with the following command:

sudo pip install sympy
sudo pip install gravipy

Then, you can import everything from both modules with the Python statements:

from sympy import *
from gravipy import *

This way, you will have access to all of the functionality from these modules within the namespace of your own Python code.

When you do work in general relativity, you often need to look at very complicated equations. The GraviPy module includes a function, called init_printing(), that sets up everything you need for pretty-printed equations. You probably will want to call this near the top of your code so that your worksheet is more human-readable.

Let's start with the SymPy Python module. SymPy is designed to give you the ability to do symbolic mathematical computations. With it, you can do things like solve algebraic expressions, rearrange and simplify equations, and even perform symbolic derivatives and integrals. I've looked at SymPy in a previous issue of LJ, so here, I just focus on some of the core parts as a reminder.

The first necessary step is to define the symbolic variables that you are going to use in your calculations. Using the symbols command, you can define them with:

x,y,z = symbols('x y z')

This way, you define the variables that define the three space dimensions in Cartesian coordinates.

If you want to have all four space-time coordinates defined in spherical coordinates, you can use this:

t, r, theta, phi = symbols('t, r, \\theta, \phi')

Remember that you are feeding a string into the symbols function. This means you need to escape any special characters to get the results you are expecting.

These commands will give you the symbolic variables that you can use in expressions. But, general relativity has a special class of variables, called coordinates, that are used to define the space-time itself. There is a class, called Coordinates, that helps define this. Using the spherical coordinates above, you can create the coordinates with the statement:

x = Coordinates('\chi', [t, r, theta, phi])

The four space-time coordinates are stored in the object x. You can access the individual coordinates by using an index. In general relativity, there are two different ways of indexing variables: covariant and contravariant indexes. To look at the element values for the covariant version, you need to use positive index values. Negative index values will give you the contravariant versions. As an example, if you wanted to get the time variable, you would use the following:

x(-1)

Right away, you should notice that this implies that GraviPy uses 1-based indexing rather than 0-based indexing.

Now that you have a set of variables to define the space-time coordinates to be used, let's move on to actual general relativity.

The core part of general relativity is the metric. The metric defines how space-time is shaped. This space-time shape is what defines the gravitational force. The usual phrase that explains what is happening is that matter tells space-time how to bend, and space-time tells matter how to move. To define the metric within GraviPy, you need to start by creating a metric tensor object. You can do so with this statement:

Metric = diag(-(1-2*M/r), 1/(1-2*M/r), r**2, r**2*sin(theta)**2)

In this example, you'll notice a new variable, named M, which represents the mass of the matter that is creating this space-time distortion. If it is an item that will remain static, you don't need to do anything extra. But, there is no reason that it should remain static. If it is something that can change, you need to use the symbols command to define it. In its most general form, the metric tensor is a 4-by-4 matrix of elements. The above example is of a diagonal metric where only the non-zero elements are along the diagonal. Once you have this tensor object, you can define the metric based on it with the statement:

g = MetricTensor('g', x, Metric)

This example gives the function the metric definition and the coordinate system to be used.

To access the various elements, you can use indices that follow the same format as above. For example, you could use the following:

g(1, 1)  ->  2M/r-1

For both vectors and tensors, you can use a special index called All that will give you every possible value for the index in question. For example, you can get the entire list of coordinates with this:

x(-All)  ->  [t r theta phi]

You can get all of the elements of the metric tensor with the statement:

g(All, All)

Now that you have a set of coordinates and a metric tensor, there are a number of standard tensors that need to be calculated to help work out what the geometry of space-time actually looks like and how things like light beams travel through this geometry. The GraviPy module includes the tensor classes Christoffel, Riemann, Ricci, Einstein and Geodesic. The Christoffel, Riemann and Ricci tensors depend only on the metric. Therefore, they all have very similar forms to create new instances and get results out of them. For example, you can define the Christoffel tensor values with this statement:

Ga = Christoffel('Ga', g)

You can get individual elements with indices, just like with the metric tensor. But, some of these tensors can have a number of uninteresting zero entries. So, you can get the non-zero values with the statement:

Ga.components

This returns a dictionary where the keys are the coordinate sets for where this particular value is located, as well as the actual non-zero value at the point.

The Einstein tensor is the one used in the actual Einstein equations for space-time, and they are a little different. In order to calculate them, you first need to calculate the Ricci tensor with this statement:

Ri = Ricci('Ri', g)

Once you have this tensor, you can calculate the Einstein tensor with this:

G = Einstein('G', Ri)

Before I leave off, I should look at two techniques that are absolutely necessary to doing general relativity. The first is index contraction. This is where you end up summing values over two of the indices in a tensor. In Python, you can do this by explicitly summing with:

Rm = Riemann('Rm', g)
ricci = sum([Rm(i, All, k, All)*g(-i, -k) for i, k in 
 ↪list(variations(range(1, 5), 2, True))], zeros(4))

These two lines are equivalent to the above single creation of the Ricci tensor. In many cases, complicated calculations are not simplified automatically. This means that you need to do this simplification explicitly with the command:

ricci.simplify()

The other important technique is to calculate geodesics, which are equations that define how light beams travel in this space-time. You need to create a new variable to handle the world line parameter for these equations:

tau = symbols('\\tau')

Now, you can calculate the geodesics with this:

w = Geodesic('w', g, tau)

Again, you can use 1-based indices to access the various non-trivial equations for your space-time of interest.

With SymPy, you can do all sorts of symbolic calculations normally reserved for programs like Maple, Maxima or Mathematica. Building on these capabilities, GraviPy lets you play with space-time and do calculations to determine curvature and gravitational effects. There is not a great deal of information available on-line explaining what GraviPy can do. Your best option is to download the source files, which include a tutorial iPython notebook, even if you install GraviPy through pip. Now you can do your gravitational research all from the comfort of Python. If you are a Python fan, this module will let you do interesting research work in your favourite language.