Tensorial calculations - Best computational Language

[Arman777]

I have given a metric in terms of the spherical coordinates $(r, \theta, \phi)$. The problem is that I need to find the Riemann/Ricci Tensor and Riccl Scalar from the given metric by using code. I wonder what type of language (such as Python, Mathematica, Matlab, etc.) I should use it to make tensorial calculations and taking derivatives easier.

 

[Ibix]

I use Maxima, which is free and multiplatform, so I have a copy on my laptop and on my phone. Take a look at my Insight article for some code to scavenge.

I believe @pervect is a Maple user. I keep meaning to see if python's sympy library has tensor/GR support, but haven't got around to it.

 

[Arman777]

I use Maxima, which is free and multiplatform, so I have a copy on my laptop and on my phone. Take a look at my Insight article for some code to scavenge.

I believe @pervect is a Maple user. I keep meaning to see if python's sympy library has tensor/GR support, but haven't got around to it.

I ll look at it. I believe that python is not a good choice for this type of problem.

 

[Ibix]

I ll look at it. I believe that python is not a good choice for this type of problem.

Base python would be a bad idea, yes. But sympy is a python symbolic maths package, so it should have the right basic tools for the job. However, I don't know if it has canned routines for generating the tensors of curved spacetime. Probably not.

Maxima does have canned calculations in the ctensor library (watch the index order though - mcs[i,j,k] is $\Gamma^k_{ij}$, and riem[i,j,k,l] is $R^l{}_{ijk}$). Maple has canned calculations too, in GRTensor, about which I know no more than the name.

 

[strangerep]

[...] I don't know if [sympy] has canned routines for generating the tensors of curved spacetime. Probably not.

I've been using sympy for quite a while now. For ordinary sympy there IS an EinsteinPy package. I don't use EinsteinPy myself, since I want to do Finsler geometry, hence I must write lower level sympy workfiles.

SymPy is best used (imho) within the jupyter-lab framework (NOT jupyter-notebook which seems not being maintained much anymore -- maybe because jupyter-lab is MUCH better).

At first, it took me quite a while to assemble all the various sub-packages, with correct versions, that are needed to make everything work properly. Once I understood how python modules fit together, and how to work with multiple workfiles in jupyter-lab, things gradually became easier.

@Arman777: What platform are you working on? I use OpenSUSE Leap. Yast2 offers jupyterlab with python3 and python3-sympy all ready for easy installation. Much of my initial difficulty was from not realizing that I should be using jupyter-lab rather than jupyter-notebook.

SymPy sometimes needs a helping hand when simplifying very complicated expressions, and it takes a while to master all the various routines that can potentially help with that. E.g., routines for splitting out numerators and denominators, simplifying each separately, and then putting them back together.

I also noticed that upgrading my (old-ish) laptop to SSD + 16GB RAM sped sympy up quite a bit.

Maple has canned calculations too, in GRTensor, about which I know no more than the name.

If one already has a Maple installation, then certainly try their GRtensor stuff (as well as whole lot of other physics-oriented packages).

I tried very hard to use Maple to help with my Finsler geometry stuff but eventually gave up. I find Maple's GUI interface more tedious than SymPy-in-Jupyter-Lab. As for computation features, my experiences were both good and not-so-good with both packages. Both Maple and SymPy have strengths and weakenesses, and navigating them in sufficient depth to assess them correctly can consume a lot of time.

Any symbolic math package requires a substantial amount of time investment to familiarize oneself with its capabilities and foibles. I'd say that time invested in SymPy (and python3 if you don't already know python) is a good investment.

[For future reference: there's also a project "symengine" which is intended to be a much faster C++-based replacement for SymPy. I get the feeling it's not yet ready (as of Mar 2021) for primetime, but I have high hopes for its future.]

 

[Ibix]

Any symbolic math package requires a substantial amount of time investment to familiarize oneself with its capabilities and foibles.

Definitely. Especially the "foibles" bit.

 

[Arman777]

Our teacher suggested Maple and Mathematica, but it seems Maple is not free, and from your suggestions, it looks Maxima/Python is a better choice.

've been using sympy for quite a while now. For ordinary sympy there IS an EinsteinPy package. I don't use EinsteinPy myself, since I want to do Finsler geometry, hence I must write lower level sympy workfiles.

Thanks, I'll look into it.

SymPy is best used (imho) within the jupyter-lab framework (NOT jupyter-notebook which seems not being maintained much anymore -- maybe because jupyter-lab is MUCH better).

Thanks for the tip as well. I was using VS Code, but for some jobs, it seems better to use Jupyter.

What platform are you working on?

Windows 10

Any symbolic math package requires a substantial amount of time investment to familiarize oneself with its capabilities and foibles.

I have never used a symbolic math package :)

I know python 3 in general and worked with available packages such as NumPy, matplotlib, scipy, etc.

The other problem is that I only have two weeks :) So I really should write the code fast and just run. So I need to have the right tools for the language and don't bother much with adjusting the tensors etc.

 

[Ibix]

Well, Maxima (or wxMaxima for a nicer interface) plus a read of my Insight article on solving the EFEs in Maxima (and Chris Hillman's older post linked therein) should get it done fairly quickly. Generating the tensors is trivial, but if they're algebraically nasty some mucking around with simplification functions may be necessary.

I'd definitely consider looking at the python packages @strangerep discussed. I suspect they're under more active development than Maxima, and python is definitely a nicer language than what is essentially a bunch of acros on top of Lisp.

 

[Arman777]

The problem is only for solving the given metric. Furthermore, the given metric equations are so simple and do not need any trick. For instance one of them is like

$$g_{\mu \nu} = \begin{bmatrix}
-1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & r^2 & 0 \\
0 & 0 & 0 & r^2sin(\theta)^2 \\
\end{bmatrix}$$

I ll try to look both and see which is nicer.

 

[PeterDonis]

one of them is like

$$g_{\mu \nu} = \begin{bmatrix}
-1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & r^2 & 0 \\
0 & 0 & 0 & r^2sin(\theta)^2 \\
\end{bmatrix}$$

You should be able to solve this one by inspection. (Hint: What does it look like in Cartesian coordinates? Second hint: Look at just the $g_{tt}$ and $g_{rr}$ metric coefficients; what do they remind you of?)

 

[Arman777]

You should be able to solve this one by inspection. (Hint: What does it look like in Cartesian coordinates? Second hint: Look at just the $g_{tt}$ and $g_{rr}$ metric coefficients; what do they remind you of?)

Oh the problem is not about solving the problem by hand but via computer. The given metric is for flat space in terms of spherical coordinates.

 

[Arman777]

I am thrilled to announce that I manage to run the code and see the output. Here are the results,

I only used 10-15 lines of code (if you include imports etc., it is about 15, otherwise it's 10). It's nice to see that the code works fine. I cannot share the code since it's homework for us, but I can share the solutions so you guys can check (well actually nothing to check since the metric corresponds to a flat universe, everything will be 0). I first started from Python with the EinsteinPy package since I know Python 3. I also used VSCode/Jupyter Notebook/Jupyter Lab, and in all of them, the code works fine. But ıt works best in Jupyter Lab/Notebook

Thank you all for your suggestions/help.

 

[vanhees71]

Hm, I'm again motivated to start learning Python...

 

[George Jones]

I have worked at three institutions over the last 15 years, and at each place Maple has been freely available to me, so I have used Maple for GR (and other stuff). At a meeting a few days ago, I was told that our site license will be maintained for one more year, but no guarantees after that. I have to get a bunch of Maple users to make a bunch of noise.

 

[Arman777]

Guys I have a problem. It seems that I need to do each step by hand and cannot use the build in functions so. I have to find the Chris. symbols by hand via loop and then same for the Riemmen tensor etc. So Is there a way to calculate the derivative of a given metric via einsteinpy or if not how can I do it...?

That is I need to find a way to calculate $\partial_{\nu}g_{\alpha \beta}$ etc.

Edit: I have found it, but it seems the code will be much longer now :)

https://www.tutorialspoint.com/sympy/sympy_derivative.htm#:~:text=we can find the differentiation,() function in SymPy package.&text=To take multiple derivatives, pass,a number after the variable.&text=It is also possible to call diff() method of an expression.

 

[Ibix]

Each element can be computed by hand. $\partial_\mu g_{\nu\rho}$ is just the derivative of the $\nu,\rho$ element of a 2d array with respect to the $\mu$th coordinate. You need to find out how to differentiate, and how to identify the symbol for the $\mu$th coordinate.

I know you aren't using Maxima, but in that language the diff(y,x) function calculates $\frac{\partial y}{\partial x}$, the ct_coords array holds the coordinate names, and the metric is in lg. So diff(lg[nu,rho],ct_coords[mu]) would do the job. You need to find how to do that in python.

 

[Arman777]

I am trying

 

[Arman777]

Each element can be computed by hand. $\partial_\mu g_{\nu\rho}$ is just the derivative of the $\nu,\rho$ element of a 2d array with respect to the $\mu$th coordinate. You need to find out how to differentiate, and how to identify the symbol for the $\mu$th coordinate.

I know you aren't using Maxima, but in that language the diff(y,x) function calculates $\frac{\partial y}{\partial x}$, the ct_coords array holds the coordinate names, and the metric is in lg. So diff(lg[nu,rho],ct_coords[mu]) would do the job. You need to find how to do that in python.

I manage to do it. The only problem is that the

 init_printing()

is not working properly for an array of the type

sympy.tensor.array.sparse_ndim_array.MutableSparseNDimArray

In this case, ,t looks like this

[[[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]], [[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, -1.0*r, 0], [0, 0, 0, -1.0*r*sin(theta)**2]], [[0, 0, 0, 0], [0, 0, 1.0/r, 0], [0, 1.0/r, 0, 0], [0, 0, 0, -1.0*sin(theta)*cos(theta)]], [[0, 0, 0, 0], [0, 0, 0, 1.0/r], [0, 0, 0, 1.0*cos(theta)/sin(theta)], [0, 1.0/r, 1.0*cos(theta)/sin(theta), 0]]]

but I want it to be like

you guys know a way ?

 

[Arman777]

I manage to do it. The only problem is that the

 init_printing()

is not working properly for an array of the type

sympy.tensor.array.sparse_ndim_array.MutableSparseNDimArray

In this case, ,t looks like this

[[[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]], [[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, -1.0*r, 0], [0, 0, 0, -1.0*r*sin(theta)**2]], [[0, 0, 0, 0], [0, 0, 1.0/r, 0], [0, 1.0/r, 0, 0], [0, 0, 0, -1.0*sin(theta)*cos(theta)]], [[0, 0, 0, 0], [0, 0, 0, 1.0/r], [0, 0, 0, 1.0*cos(theta)/sin(theta)], [0, 1.0/r, 1.0*cos(theta)/sin(theta), 0]]]

but I want it to be like

View attachment 280871

you guys know a way ?

I have found it!

 

[PeterDonis]

I have found it!

So what is it?

 

[PeterDonis]

I cannot share it here

Yes, you can. What's more, you should. Otherwise you are depriving other readers of this thread of the answer to your question. If you want to post questions here, that's part of the deal.

 

[Ibix]

You don't need to post your complete solution (I presume you don't want to as it's homework). But whatever fix you did to init_printing to make it work with the sparse array type you were using would be useful to know.

 

[strangerep]

Regarding init_printing(), I always put something like this at the start of my SymPy files:

import sympy
from sympy import *

sympy.init_printing()

It helps to learn a bit about how to define and use your own Python modules. It's pretty easy once you grok the pattern. E.g., create a MyModule.py file. Then, in your ipynb (file -- i.e., the file that jupyter-lab will actually use), put:

import MyModule
from MyModule import *

which imports all the names, functions, etc, from MyModule. Without the "from" statement, you'd need to say MyModule.somefunc(...), or whatever, instead of just somefunc(...).

Jupyter-lab seems reasonably smart in that it will cache a pre-munched version of MyModule.py in a "__pycache__" subdirectory, which helps to speed up subsequent re-runs.

Regarding use of EinsteinPy:

https://docs.einsteinpy.org/en/latest/examples/Symbolic%20Module%20-%20Einstein%20Tensor.html

If you poke around in that examples area, it also shows use of RicciTensor, etc, etc.

(Warning: some of the EinsteinPy tutorials were originally written to run in jupyter-notebook. Now they give me "Error 404" errors. Running EinsteinPy locally in jupyter-lab would probably work much better -- but I haven't actually tried it.)

Regarding direct use of the tensor classes in SymPy, I could never get them to work properly. I ended up using implicit arrays of arrays, etc. It also helps to learn how to define and use basic functions in Python. E.g.,

# Given an N^X_x, and velocity 4-vector, compute the corresponding B^X_{xy}
# and return it as a matrix.
#
def mmBXxx(NXx, uv):
    result_BXxx = [
                                 [[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]],
                                 [[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]],
                                 [[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]],
                                 [[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]]
                           ]
for X in range(4):
    for x1 in range(4):
        for x2 in range(4):
            result_BXxx[X][x1][x2] = NXx[X][x1].diff(uv[x2])

return result_BXxx

(NXx is 2nd-rank, BXxx is 3rd-rank, uv is a 1D array of variables.)

Probably, I should study the internals of the EinsteinPy module properly to see how the experts do it.

[Edit: There's a lot of code in the EinsteinPy module, and it drags in a LOT of other stuff because it's trying to be all things to all,... er,... relativists. I think I'll stick with writing my own lower level, but much simpler, code for now.]

 

[strangerep]

I have worked at three institutions over the last 15 years, and at each place Maple has been freely available to me, so I have used Maple for GR (and other stuff). At a meeting a few days ago, I was told that our site license will be maintained for one more year, but no guarantees after that. I have to get a bunch of Maple users to make a bunch of noise.

Try SymPy and EinsteinPy first.

 

[vanhees71]

I have worked at three institutions over the last 15 years, and at each place Maple has been freely available to me, so I have used Maple for GR (and other stuff). At a meeting a few days ago, I was told that our site license will be maintained for one more year, but no guarantees after that. I have to get a bunch of Maple users to make a bunch of noise.

The same with me. Since my diploma thesis times at all institutions we had Mathematica, and that's why I'm used to Mathematica for symbolic computer algebra. I think it's pretty common that one uses the tools one has access to and then stays with them, because one is used to them.

The problem with commercial software is of course that particularly for institutions like universities or other research institutions the licences are particularly expensive. So one has regularly debates about whether to prolong the license agreements or even buy more :-(.

 

[Arman777]

Yes, you can. What's more, you should. Otherwise you are depriving other readers of this thread of the answer to your question. If you want to post questions here, that's part of the deal.

Sorry guys. Okay so this works for me

I just turned it into an array by using

.tolist()

Such as

Christoffel_symbol.tolist()
Christoffel_symbol

If you poke around in that examples area, it also shows use of RicciTensor, etc, etc.

I did. That site helped me. I cannot use Einsteinpy since that's not what our teacher wants :p. But I have calculated all of the tensors (Riemann, Ricci) and Ricci scalar without using Einsteinpy. All you guys need is sympy and NumPy. It's not hard to code, since there are only some "for" loops. But it takes some time to write.

Regarding direct use of the tensor classes in SymPy, I could never get them to work properly. I ended up using implicit arrays of arrays, etc. It also helps to learn how to define and use basic functions in Python. E.g.,

Instead of tensors, I have used arrays such as

xtensor = MutableSparseNDimArray(zeros((4,)*2))

By doing this, you can allow for updating your matrix and changing the elements. Otherwise, the code was giving me errors.
Maybe after the homework deadline passes, I can share the code.

 

[MathematicalPhysicist]

I have worked at three institutions over the last 15 years, and at each place Maple has been freely available to me, so I have used Maple for GR (and other stuff). At a meeting a few days ago, I was told that our site license will be maintained for one more year, but no guarantees after that. I have to get a bunch of Maple users to make a bunch of noise.

There's nothing like the holy-M trinity (Matlab+Maple+Mathematica); pricey I know.

 

[Arman777]

Here you can see the full solution in a code or as a pdf. I appriecete any suggestion or improvements. I am also thinking turning this into a class and make it something like Einsteinpy but ofc maybe I can do that in summer.

https://github.com/seVenVo1d/General-Relativity-Tensor-Calculations

 

[MathematicalPhysicist]

CHEERS MATE

 

[strangerep]

Here you can see the full solution in a code or as a pdf. I appriecete any suggestion or improvements.

Try a general spherically symmetric (but still diagonal) metric -- e.g., suppose you didn't know the Schwarzschild metric and are trying to find it. Can you make the Einstein equations reduce down to a set of differential equations that's fairly easy to solve?

If so, then improve your code to handle non-diagonal metrics (you'll need to use a matrix inverse function somehow to get the inverse metric). Can you set up the Einstein equations for a general axi-symmetric metric, and simplify it down to a set of DEs that can be solved to find the Kerr metric?

I am also thinking turning this into a class and make it something like Einsteinpy but ofc maybe I can do that in summer.

https://github.com/seVenVo1d/General-Relativity-Tensor-Calculations

EinsteinPy tries to be all things to all people, but fails to be good at some things which are relatively difficult (like simplifying the Einstein equations to get a mangeable set of DEs). I suggest keeping different capabilities in distinct packages (maybe using a common set of library code if applicable). I.e., avoid the mega-tar-baby approach to coding.

 

[Arman777]

.g., suppose you didn't know the Schwarzschild metric and are trying to find it. Can you make the Einstein equations reduce down to a set of differential equations that's fairly easy to solve?

Well, for simple equations, maybe the computer can solve some problems...or I can pre-upload some solutions and metrics, etc.

If so, then improve your code to handle non-diagonal metrics (you'll need to use a matrix inverse function somehow to get the inverse metric). Can you set up the Einstein equations for a general axi-symmetric metric, and simplify it down to a set of DEs that can be solved to find the Kerr metric?

like simplifying the Einstein equations to get a mangeable set of DEs)

I have never tried that thing, but maybe it's due to problems of the sympy. Sadly it seems that the sympy does not have enough power to simplify expressions as Mathematica. So even simple terms remain complicated due to this issue, which I have encountered.

I updated my code (https://github.com/seVenVo1d/General-Relativity-Tensor-Calculations/blob/main/GR_tensorial_cal_v2.ipynb) and add Weyl Tensor Einstein Tensor, and also two forms of the Riemannian Tensor ($R^{\alpha}_{\beta \gamma \nu}, R_{\alpha \beta \gamma \nu}$).

I was thinking of including index lowering and raising function and maybe improving the style of the code.

I suggest keeping different capabilities in distinct packages (maybe using a common set of library code if applicable). I.e., avoid the mega-tar-baby approach to coding.

I definately need to do that. Either I am going to create a class which I need to do at this point. In that case it would be easy to treat matrices as objects but I am kind of new to creating classes s that would take some time. After creating the classes it would be easy to work with the code I guess. Even in this case I could create another .py file and just use the functions but I'll keep this way for now.

 

[strangerep]

...or I can pre-upload some solutions and metrics, etc.

I don't see much point in that.

I have never tried that thing, but maybe it's due to problems of the sympy. Sadly it seems that the sympy does not have enough power to simplify expressions as Mathematica. So even simple terms remain complicated due to this issue, which I have encountered.

It might be interesting to try a particular simplification task in SymPy, Maple and Mathematica and see which one can progress furthest without human assistance. And then,... with a little bit of human assistance...

I was thinking of including index lowering and raising function and maybe improving the style of the code.

I think you need to clarify in your own mind what real-world problem(s) you are really trying to solve by creating this software. If you produce some nicely structured, well-written code that's not actually useful in the real world, then what's the point?

 

[Arman777]

I think you need to clarify in your own mind what real-world problem(s) you are really trying to solve by creating this software. If you produce some nicely structured, well-written code that's not actually useful in the real world, then what's the point?

I thought a user would enter or type the metric components, and then it will have access to Riemann Tensor, Ricci Tensor, Einstein Tensor, etc. Then, he/she can raise or lower the indices of these tensors at any level. This was my purpose. Solutions of the GR equations only can be possible after these steps, I guess. To do these things, I need to write a class and separate the function in one .py file and maybe try to write a user-friendly interface instead of typing the code by hand.

It might be interesting to try a particular simplification task in SymPy, Maple and Mathematica and see which one can progress furthest without human assistance. And then,... with a little bit of human assistance...

I believe it's Mathematica...I just need to add some of the simplifaction functions so that the code gets simplified. And sometimes simplification makes things worse in sympy.

I don't see much point in that.

Well, maybe not much but just a little ... I am doing/want to do things done by Einsteinpy, but I am focusing on the curvature tensors rather than solving tensor equations.

 

[Arman777]

I manage to turn it into a class. I have also put a detalied description about how to run the code