Christoffel symbols of the first kind ill calculated - Calculation of the Christoffel symbols of the second kind returns an error

[FreddyBaudine]

Dear All,

Galgebra does not seem to obtain the correct Christoffel symbols of the first kind; at least the values printed out indicate a skew symmetry, which is strange.

Starting with the metric tensor relating to a two dimensional manifold, the unit sphere, viz., galgebra yields the following components for the Christoffel symbols of the first kind, when, in fact, the correct values are As is well known, these symbols are symmetrical in the indices alpha and beta.

As for the Christoffel symbols of the second kind they are not calculated, an error is returned instead: “ C:\Users\Freddy\Desktop\Issues\Christoffel>C:\Anaconda3\python "Christoffel_symbols_checking.py" Traceback (most recent call last): File "Christoffel_symbols_checking.py", line 47, in Cf2 = sp2.Christoffel_symbols(mode=2) File "C:\Anaconda3\lib\site-packages\galgebra\metric.py", line 490, in Christoffel_symbols Gamma_ijk += Gijl * self.g(l,) TypeError: 'MutableDenseMatrix' object is not callable

C:\Users\Freddy\Desktop\Issues\Christoffel>pause Press any key to continue . . . “ It is not at all clear how the differential operators can be properly computed (see zip file enclosed) when it is explicitly stated in Section 2.3.1 of Galgebra: a Geometric Algebra Module for Sympy (October 8, 2016) that the Christoffel symbols of the first kind are made use of to perform a number of calculations, notably to obtain the derivatives of the basis vectors (footnote on page 17). Unless the Christoffel symbols are not actually used at all! This would explain that, but I can be totally wrong on this.

Look forward to receiving any views on all this. Kind regards, FB freddy.baudine@outlook.com Christoffel.zip

[utensil]

It is not at all clear how the differential operators can be properly computed (see zip file enclosed) when it is explicitly stated in Section 2.3.1 of Galgebra: a Geometric Algebra Module for Sympy (October 8, 2016) that the Christoffel symbols of the first kind are made use of to perform a number of calculations, notably to obtain the derivatives of the basis vectors (footnote on page 17). Unless the Christoffel symbols are not actually used at all! This would explain that, but I can be totally wrong on this.

Historically, derivatives_of_basis() is implemented before Christoffel_symbols() and the latter is currently not covered by any examples or tests, see https://codecov.io/gh/pygae/galgebra/src/15-print-pow/galgebra/metric.py#L426 .

Christoffel_symbols() was implemented in commit https://github.com/rschwiebert/galgebra/commit/e9177cc128db2eb17a1fe79fb7cdb2d6285c9cf3 and not changed ever since.

And there's duplicate code in derivatives_of_basis() for calculating Christoffel symbols of the first kind. The code is almost identical and straight-forward, I'm still investigating what could have gone wrong.

[utensil]

@FreddyBaudine According to the implementation, indices need to be interpreted as Christoffel_symbols()[i, j, k] = \Gamma_{kij}, note the order of k here. In this interpretation, the result is correct and identical to the result given by Maple.

And I have confirmed that the derivatives_of_basis() has used such semantics so it's also correct.

What's wrong here, is that Christoffel_symbols()[i, j, k] didn't have a proper latex printer to reflect such semantics.

[utensil]

On the other hand the Christoffel symbols of the second kind are not implemented correctly, working on a fix.

[utensil]

@FreddyBaudine I've just recreate the same form as your Maple file by commit 303acab , it's basically just:

# Add these imports at the beginning of the Python script
from sympy import Array, permutedims

# Replace with these lines for calculating Cf1
Cf1 = sp2.Christoffel_symbols(mode=1)
Cf1 = permutedims(Array(Cf1), (2, 0, 1))
print(r'\text{Christoffel symbols of the first kind: }')
print(r'\Gamma_{1, \alpha, \beta} = ', latex(Cf1[0, :, :]), r'\quad', r'\Gamma_{2, \alpha, \beta} = ', latex(Cf1[1, :, :]))

The result looks like this:

[FreddyBaudine]

Dear Utensil, Many thanks for all this. Your doing a fantastic job. I tried your solution with Python 2.7, 3.4 and 3.7 together with version 1.3 of Sympy. Everything turned out alright. Kind regards, FB

[utensil]

Dear Freddy,

I'm glad it checks out and it's a pleasure to be helpful to your work. Also thanks for your verification on multiple platforms which are not verified on CI yet (only 2.7 and 3.6 for now).

P.S. Did you notice my invite for you to join our discussion on https://clifforddev.slack.com/ ? I sent the invitation to your email address mentioned in the main thread of this issue.

Regards, Utensil

Slack

[utensil]

Dear @FreddyBaudine ,

The calculation of Christoffel symbols of the second kind has been implemented in https://github.com/pygae/galgebra/blob/15-print-pow/galgebra/metric.py#L467 .

And an equivalent result of your script has been implemented in https://github.com/pygae/galgebra/blob/15-print-pow/examples/LaTeX/christoffel_symbols.py#L48 and it looks like this:

It's slightly different from your Maple file, because it simplified cos(theta)/sin(theta) to 1/tan(theta).

In the process, I've refactored it quite a bit, so now derivatives_of_basis() actually uses Christoffel_symbols() for calculation instead of duplicated but slightly different code.

Also I would like the code to look like the original math formulas because reading the original code style can't help me eyeball the correctness of the implementation which I believe it's much improved now. Please take a look, comments welcome!

Regards, Utensil

[FreddyBaudine]

Dear Utensil, As before, I tried your solution with Python 2.7, 3.4 and 3.7 together with version 1.3 of Sympy. Everything turned out alright. Kind regards, Freddy

[utensil]

Dear @FreddyBaudine,

Thanks for testing the solutions. Please note that galgebra currently have issues with SymPy 1.4, see #31 .

Have you verified the calculation of Christoffel symbols (both for the first and second kind) for other metrics? Although I'm sure the implementation is correct, but it's better to know that it works for more scenarios.

There's one more thing I need to ask you: During the fix of this issue, I added your code to the test cases and examples of glgebra, it's fine for an unmerged PR #17, but it might have licensing issue if they become part of galgebra. What do you think? Do we have your permission to add them into galgebra test cases and examples?

[FreddyBaudine]

Dear Utensil,

I checked your corrections with higher dimension metrics and, unfortunately, the Christoffel symbols of the second kind are not computed correctly in some cases. As far as I can judge, this mishap happens only if the matrix of the metric possesses off-diagonal elements. That’s why I tried two such metrics, viz., the rotating frame metric and the Kerr-Debney metric for they are not too complicated. I will test the Kerr metric eventually.

It appears that, in all cases envisaged, the Christoffel symbols of the first kind are computed correctly.

I did all this making use of Python 3.7 and Sympy version 1.3

Kind regards, Freddy PS: A priori, no license issue as this work is strictly personal. Hence, you can add my code in galgebra test cases and examples. 37.zip

[utensil]

Dear @FreddyBaudine,

Great catch!

Currently I have no idea where went wrong, because both the metric inverse code and the Christoffel symbols of the second kind code seem ok to me. Still investigating.

I'm cross-checking the implemtation and test cases in galgebra and sympy.diffgeom.

A priori, no license issue as this work is strictly personal. Hence, you can add my code in galgebra test cases and examples.

That would be greate, thanks a lot!

Regards,

Utensil Song

[utensil]

Also I intend to also verify the implementation against some cases in https://github.com/spacetimeengineer/spacetimeengine/blob/master/spacetimeengine/src/solutions.py .

GitHub

spacetimeengineer/spacetimeengine

A Python utility for analyzing a given solution to the Einstein's field equations. Built on Sympy. - spacetimeengineer/spacetimeengine

[utensil]

I've checked the Kerr-Debney metric in your attached script, it turns out the the Christoffel symbols of the second kind is correctly calculated. They don't look identical to the Maple one because of two reasons:

First, we can't use 1/2 to express half in SymPy but should use Rational(1, 2) instead to avoid being treated as 0.5 and also the bug that galgebra will translate it to something like 0 \cdot 5 for all decimal numbers.

You can do

from galgebra.mv import ONE, ZERO, HALF

and write the metric as:

g = [[0, 0, -exp(-z), 0],
     [0, HALF*u**2*exp(4*z), 0, 0],
     [-exp(-z), 0, 12*exp(-2*z), u*exp(-z)],
     [0, 0, u*exp(-z), HALF*u**2]]

Now it will be much easier to see the next reason.

Secondly, if we write:

Cf2 = Cf2.normalized()

Then the result is identical to the Maple one:

[FreddyBaudine]

Dear Utensil,

Very many thanks for your input. I tried to reproduce your solution but I got the following error with Cf2.normalized() :

“ C:\Users\Freddy\Desktop\Issues\Christoffel\Higher dim\37>C:\Python37\python "Christoffel_kerr_debney_metric37.py" Traceback (most recent call last): File "Christoffel_kerr_debney_metric37.py", line 48, in Cf2 = Cf2.normalized() AttributeError: 'list' object has no attribute 'normalized'

C:\Users\Freddy\Desktop\Issues\Christoffel\Higher dim\37>pause Press any key to continue . . . “ How come I got this error when you did not? For your convenience I have enclosed the modification in a zip file as I may have made some mistake. If this is the case, will you be so kind as to forgive me to have wasted your time. You will have guessed that I am not quite competent with Python. I am much more familiar with Maple, Mathematica and a little bit of maxima.

Thanks again to take the time to answer my queries.

Kind regards, FB Christoffel_kerr_debney_metric37.zip

[utensil]

@FreddyBaudine Oh I forgot to mention that normalized() should be called after permutedims() call because after that it will be an array and it will have such method.

_{Sent with GitHawk}

[FreddyBaudine]

Dear Utensil,

Many thanks for your quick reply. I placed the instruction Cf2 = Cf2.normalized() after the instruction Cf2 = permutedims(Array(Cf2), (2, 0, 1)) as you suggested. Now I get the following error :

“ C:\Users\Freddy\Desktop\Issues\Christoffel\Higher dim\37>C:\Python37\python "Christoffel_kerr_debney_metric37.py" Traceback (most recent call last): File "Christoffel_kerr_debney_metric37.py", line 49, in Cf2 = Cf2.normalized() AttributeError: 'ImmutableDenseNDimArray' object has no attribute 'normalized'

C:\Users\Freddy\Desktop\Issues\Christoffel\Higher dim\37>pause Press any key to continue . . . “

I trust you will find what’s wrong.

Kind regards, FB Christoffel_kerr_debney_metric37_Corrected.zip

[utensil]

Sorry, my mistake.

At first, I noticed galgebra's output is different by a factor of * (-4 * exp(-2 * z) * u**(-4)), so I used

Cf2 = Cf2 * (-4 * exp(-2 * z) * u**(-4))

to achieve the expected result. But I think this is caused by lack of normalization so I added the normalized() call but it's actually incorrect. And unfortunately I ran the script by repeating the last command which actually didn't run the script and only open the last pdf file generated by Cf2 = Cf2 * (-4 * exp(-2 * z) * u**(-4)).

Sorry for the inconvenience caused and I was on other computer and didn't have access to the script and gave the wrong guide based on impression.

So, I need to find the actual fix.

P.S. The factor is actually det(g)^-1.

[utensil]

I believe we should use norm=True when initializing the GA object, but it will raise an error:

  File "~/Christoffel_kerr_debney_metric37.py", line 23, in <module>
    g4 = Ga('eu ex ey ez', g = g, coords = g4coords, norm=True) # Create g4
  File "~\projects\pygae-galgebra\galgebra\ga.py", line 325, in __init__
    metric.Metric.__init__(self, bases, **kwargs)
  File "~\projects\pygae-galgebra\galgebra\metric.py", line 700, in __init__
    self.e_norm.append(square_root_of_expr(self.g[i, i]))
  File "~\projects\pygae-galgebra\galgebra\metric.py", line 108, in square_root_of_expr
    (coef, pow_lst) = sqf_list(expr)
  File "~\Miniconda3\envs\py3\lib\site-packages\sympy\polys\polytools.py", line 6216, in sqf_list
    return _generic_factor_list(f, gens, args, method='sqf')
  File "~\Miniconda3\envs\py3\lib\site-packages\sympy\polys\polytools.py", line 5988, in _generic_factor_list
    raise PolynomialError("a polynomial expected, got %s" % expr)
sympy.polys.polyerrors.PolynomialError: a polynomial expected, got 12*exp(-2*z)

I read the code a little bit, it seems to be a bug in square_root_of_expr and it seems that it will take quite some time to fix it.

Please simply times det(g)^-1 like this for now:

Cf2 = g4.Christoffel_symbols(mode=2)
Cf2 = permutedims(Array(Cf2), (2, 0, 1))
Cf2 = Cf2 * (g4.g.det() ** -1)

[FreddyBaudine]

Dear Utensil,

Thank you for the workaround.

I trust you won’t be surprised of the following. As the Christoffel symbols of the second kind are calculated correctly when the matrix of the metric is diagonal, the workaround should not work in that situation. Sure enough, this is the case; all the components get divided by the determinant. I’ve tested this with the Schwarzschild metric.

Kind regards, FB

[utensil]

The direct cause is found. The Christoffel symbols of the second kind relies only on the Christoffel symbols of the first kind and the inverse of metric. The former is correctly calculated, so I initially suspected the inverse of metric, but when looking at your pdf, the inverse of metric seems to be correct. Unfortunately, it's correct simply because yours is calculated by simply g.inv() and galgebra did something completely different so the inverse of metric used was incorrectly calculated by exactly det(g)^-1.

At first I thought the problem is in Mv.inverse_metric() but it looks logical:

    def inverse_metric(self):

        if self.g_inv is not None:
            return

        if self.is_ortho:  # Orthogonal metric
            self.g_inv = eye(self.n)
            for i in range(self.n):
                self.g_inv[i,i] = S(1)/self.g(i,i)
        else:
            if self.gsym is None:
                self.g_inv = simplify(self.g.inv())
            else:
                self.detg = Function('|' +self.gsym +'|',real=True)(*self.coords)
                self.g_adj = simplify(self.g.adjugate())
                self.g_inv = self.g_adj/self.detg
        return

And I checked the conditions, it should execute the line self.g_inv = simplify(self.g.inv()), but simplify(self.g.inv()) actually gives the correct inverse.

After some debugging and code searching, I finally found the inverse is actually calculated by Ga.build_reciprocal_basis():

        g_inv = eye(self.n)
        self.dot_mode = '|'

        # Calculate inverse of metric tensor, g^{ij}

        for i in self.n_range:
            rx_i = self.r_symbols[i]
            for j in self.n_range:
                rx_j = self.r_symbols[j]
                if j >= i:
                    g_inv[i, j] = self.dot(self.r_basis_dict[rx_i], self.r_basis_dict[rx_j])
                    if not self.is_ortho:
                        g_inv[i, j] /= self.e_sq
                else:
                    g_inv[i, j] = g_inv[j, i]

        self.g_inv = g_inv

The self.e_sq in the executed line g_inv[i, j] /= self.e_sq is just (eu^ex^ey^ez) * (eu^ex^ey^ez) and it's equal to det(g). So it should have already been divided once, why isn't it enough?

Still investigating. The GA way of doing this is quite different from ordinary GR textbooks.

[utensil]

Dear @FreddyBaudine ,

I finally figured out that it's caused by #32, please see https://nbviewer.jupyter.org/github/utensil/julia-playground/blob/master/py/galgebra_reciprocal_basis.ipynb for step-by-step investigation.

Simply put, the inverse metric calculation in Ga is wrong for non-orthogonal basis.

This is because the calculation of reciprocal basis should use E_n^1 but choose to use E_n instead so all calculation that involves reciprocal basis need to remember to divide the factor E_n^2. In the case of inverse metric, it's the dot product of two reciprocal basis vector, so the factor should be E_n^4 but implemented as E_n^2 possibly due to oversight.

So the line g_inv[i, j] /= self.e_sq should be g_inv[i, j] /= self.e_sq**2, see 08f33f9 . (The commit passes all test cases but one, I'll fix the failed case next week).

Regards,

Utensil Song

[FreddyBaudine]

Dear Utensil,

Thank you for your involvement in all this, this is most appreciated.

I downloaded your corrected version of ga.py and run my little codes all over again, this times round obtaining the correct results. While I was casting an eye on the Jupyter nbviewer, I noticed that you designed to verify the logic step by step using “Rotating frame spacetime manifold (Kerr-Debney Metric)”. I am afraid I made a mistake in as much as the Kerr-Debney metric has got nothing to do with the rotating frame metric. Nothing very unfortunate, I simply forgot to delete the two print instructions in the file I originally posted, viz., “Christoffel_kerr_debney_metric37.py”. Sorry about that, corrected file enclosed.

In your message, you’re explaining that the commit passed all tests but one. May I ask you which one that is?

Kind regards, FB ga_corrected.zip

[utensil]

Dear @FreddyBaudine,

The Kerr-Debney Metric typo is fixed, thanks for the heads up.

In your message, you’re explaining that the commit passed all tests but one. May I ask you which one that is?

examples/LaTeX/lin_tran_check.py failed and it's now fixed (in the sense of the mismatched output, not-simplified g_inv and poorly printed matrices) at https://github.com/pygae/galgebra/pull/17/commits/84c06ed24676feb0a9db6abfaa6535d17d2bcf66 . But it always had issues: it's not properly simplified using GA identities so the formulas are very loooooong.

I also tried some other metrics in https://github.com/spacetimeengineer/spacetimeengine/blob/master/spacetimeengine/src/solutions.py at https://nbviewer.jupyter.org/github/utensil/julia-playground/blob/master/py/galgebra_gr_metrics.ipynb , they seem to have correct results. ( BTW the project is incomplete and has quite some typos.)

I have tried the Kerr metric, but GAlgebra has serious performance issue and it takes more than 20 minutes to initialize GA (I finally gave up trying for now until I get to fix #23 which I have no clue yet besides delay the simplification to the printing stage ).

In the process I also find some annoying bugs in how GAlgebra handles printing, particularly for lists, Arrays and Matrices, I'll try to fix them some time later because they likely require substantial structure refactor.

P.S. May I ask what's the background of your work on metrics with Maple and GAlgebra? More precisely, I wonder what more do you need from GAlgebra? Because in the process of fixing these issues, I realize that GAlgebra is far from complete for treating GR problems:

• it doesn't calculate Ricci tensor or Riemann curvature tensor etc.
• it uses semi-traditional semi-GA method to do certain calculations
• it doesn't make use of full power of GA method (i.e. the various identities to facilitate simplifying or calculation) And I would like to make GAlgebra (by modifying it or implement something based on it) more suitable for your use case (if feasible).

Regards,

Utensil Song

GitHub

spacetimeengineer/spacetimeengine

A Python utility for analyzing a given solution to the Einstein's field equations. Built on Sympy. - spacetimeengineer/spacetimeengine

Jupyter Notebook Viewer

[FreddyBaudine]

Dear Utensil,

Sorry for my late reply but I had been taken ill. In one of your messages you had requested that I verified the calculation of the Christoffel symbols for other metrics, possibly on higher dimensional manifolds. I simply found it easy to use some GR metrics that were readily accessible.

I understand that the task ahead is huge if Galgebra is ever to be used to solve advanced problems.

Kind regards, Freddy

[utensil]

Dear FreddyBaudine,

Sorry for my late reply, I was under the wrong impression that I've replied.

I had been taken ill

I'm sorry to hear that. Have you become better now?

Although I haven't been working on galgebra directly for a while, mostly because of busy work, I'm actively learning and preparing for further work on galgebra. Also, I've written a Julia wrapper of galgebra at https://github.com/pygae/GAlgebra.jl with fancier syntax using Unicode symbols as infix, prefix and postfix operators. Just mention it in case you like Julia too.

We had a lot of discussions in https://clifforddev.slack.com/ , and it would be great to hear from you on the discussions.

Regards,

Utensil Song

GitHub

pygae/GAlgebra.jl

Julia interface to GAlgebra via PyCall. Contribute to pygae/GAlgebra.jl development by creating an account on GitHub.

Slack