General Relativity in Python
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
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:
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
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
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
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:
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
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
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:
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.
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
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
tensor. In many cases, complicated calculations are not simplified
This means that you need to do this simplification explicitly
with the command:
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.