Kleiner SymPy-basierter Christoffel-Rechner für Python


Die Relativitätstheorie basiert auf dem Relativitätsprinzip von Galileo. Hier ist Raum, Fragen zur Relativität zur diskutieren.

langlebiges Kaon

Beiträge: 444

Registriert: Mi 26. Okt 2011, 15:34

Beitrag Mo 10. Feb 2014, 18:59

Kleiner SymPy-basierter Christoffel-Rechner für Python

Ich war's heute leid, mal wieder Christoffel-Symbole zu Fuss auszurechnen; und Maxima, Maple et al. haben dafür m.E. zuviel Overhead.

Das ist ein v0.1 quick and dirty hack von heute; ohne Gewähr jedweder Art; Doku im Code bitte beachten!

  Code:
replaced with version 0.1.1 below


  Code:
#!/usr/bin/python2.7
# -*- coding: utf-8 -*-
#
#  christoffel_calculator.py

#  Copyright 2014 Solkar (Sk)
#  <http://www.quantenforum.de/memberlist.php?mode=viewprofile&u=482>
#
#  version history:
#       0.1 (buggy and widely untested) - Feb 10th 2014 (Sk)   
#            @BUG yet just one possible dump per run
#                 (constraint realized by main() anyway)
#                 likely this a side effect related to sympy symbol allocation
#            @TODO implement simpification of generated symbolics if possible
#            @TODO design a class for encapsulation of the yet globals
#            @TODO add means to attach symbolic indices t,r.. instead of 0,1.. 
#                  to the Γ symbols 
#       0.1.1 (still buggy and widely untested) - Feb 11th 2014 (Sk)   
#            Fixed a typo in FLRW metric, improved Unicode output
#            BUGs and TODOs as before
#
###############################################################################

#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 of the License, or
#  (at your option) any later version.

#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.

#  You should have received a copy of the GNU General Public License
#  along with this program; if not, write to the Free Software
#  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
#  MA 02110-1301, USA.

###############################################################################
"""
    sympy¹-based snippets for generating symbolic representations of
    Christoffel symbols

    GR-related conventions:
    c := 1, t leading, index numbering starts at 0, +2 signature

    ¹http://sympy.org/en/index.html

    Bibliography:
    [MTW73] Misner, C.; Thorne, K. & Wheeler, J. Gravitation.
            W. H. Freeman, 1973.
"""

import sys       # for argv, stderr 
import itertools # for index tuple generation
import codecs
UTF8Writer = codecs.getwriter('utf8')
sys.stdout = UTF8Writer(sys.stdout)

import sympy
import sympy.diffgeom as dg
import sympy.printing.pretty as pp

DIM = 4
CHRISTOFFEL_INDEX_TUPLES = itertools.product(range(DIM), repeat = 3)

LATEX_ENCODER = lambda t,C: '\Gamma_{' + str(t[1]) + str(t[2]) + '}^{' + str(t[0]) + '} = ' + sympy.latex(C)
#UTF8_ENCODER = lambda t,C: 'Γ^' + str(t[0]) + '_' + str(t[1]) + str(t[2]) + " = " + repr(C)
UTF8_ENCODER = lambda t,C: pp(sympy.Symbol(u'\u0393' + "^" + str(t[0]) + '_' + str(t[1]) + str(t[2])),wrap_line = False) + unicode(" =\n") + pp(C)

# + " = "

ENCODER_LABELS=['utf-8','latex']
ENCODERS = {s:f for s,f in zip(ENCODER_LABELS, [UTF8_ENCODER, LATEX_ENCODER])}

def dump_christoffels(christoffels, encoder):
    '''
        delegates output of non-vanishing Γ^i_jk to <encoder>
    '''
    for (i,j,k) in CHRISTOFFEL_INDEX_TUPLES:
        C = christoffels[i][j][k]
        if C != 0:
            print(encoder((i,j,k),C))

   


M = dg.Manifold('M', DIM)
SSP_CS = dg.CoordSystem('Schwarzschild_polar', dg.Patch('P', M),
            ['t', 'r', 'theta', 'phi'])

#
# it's not exactly nice polluting the global namespace with lowercase
# symbols which are also CS-dependent, but python does neither provide
# file-local mods like Perl nor a "using" directive like C++11
#
# @TODO Design a singleton class for encapsulation
#
t, r, theta, phi = SSP_CS.coord_functions() 
dt2, dr2, dtheta2, dphi2 = [dg.TensorProduct(dx,dx)
                                    for dx in SSP_CS.base_oneforms()] 
#

def ds2_FLRW_SSP():
    """
        returns symbolic representation
        of FLRW line element (see e.g. [MTW73] eqns(27.24,27.27))
        in Schwarzschild polar coordinates
    """
    k = sympy.symbols('k')
    dsigma2 = ((1/(1-(k*r**2))) * dr2 + r**2 * (dtheta2 + (sympy.sin(theta))**2 * dphi2))
    a = sympy.symbols('a',cls=sympy.Function)
    ds2 = -dt2 + a(t)**2 * dsigma2
    return ds2

def ds2_SS_SSP(use_Rs = True):
    '''
        returns symbolic representation
        of Schwarzschild line element
        in Schwarzschild polar coordinates

        use_rs determines whether "2*M" or "R_s" is used for the symbolics
    '''
    ds2 = None
    f = None
    if use_Rs:
        R_s = sympy.symbols('R_s')
        f = lambda r : 1-R_s/r
    else:
        m = sympy.symbols('m')
        f = lambda r : 1 - 2*m/r
    ds2 = -f(r)*dt2 + 1/f(r) * dr2 + r**2 * (dtheta2 + (sympy.sin(theta))**2 * dphi2)
    return ds2

DS_LABELS = ('FLRW', 'Schwarzschild')
DS_MAKERS = {s:f for s,f in zip(DS_LABELS, [ds2_FLRW_SSP,ds2_SS_SSP])}


def usage():
    s = "usage: python2.7 christoffel_calculator.py ("
    s += "|".join(DS_LABELS)
    s += ") ("
    s += "|".join(ENCODER_LABELS)
    s += ")\ne.g.:"
    s += "\n\tpython2.7 christoffel_calculator.py Schwarzschild utf-8"
    return s


def main():
    """
        @TODO check why just single invocation (as enforced here) runs are
              possible;
              invoking dump_christoffels(..) more than once with
              e.g. a different metric would quietly fail on the #>1 run
        @BUG  dto.
    """
    try:
        # this is MEANT to throw on invalid invocation following
        # "EAFP" paradigma
        # "Easier to Ask for Forgiveness than Permission"
        # see http://docs.python.org/2/glossary.html
        metric_name = sys.argv[1]
        encoding = sys.argv[2]
        ds2 = DS_MAKERS[metric_name]()
        encoder = ENCODERS[encoding]           

        dump_christoffels(dg.metric_to_Christoffel_2nd(ds2), encoder)
       
    except (IndexError, KeyError) as e:
        sys.stderr.write(usage() + "\n")

if __name__ == '__main__':
    main()


EDIT: Feb 11th 2014 00:45 CET applied a patch
Old MacDonald had a form
$e_i\,\wedge\,e_i = 0$
(N.N. on Physics Forums)

Anti-Neutron

Beiträge: 921

Registriert: So 16. Jan 2011, 18:11

Beitrag Di 2. Dez 2014, 09:27

Hallo Solkar,

ich habe das Skript mal etwas erweitert, so dass man mit dem Befehl "python christoffel_calculator.py deSitter" auch die Christofffel-Symbole für diese Metrik: http://de.wikipedia.org/wiki/De-Sitter- ... oordinaten bekommt

  Code:
new version 0.1.2 below


  Code:
#!/usr/bin/python2.7
# -*- coding: utf-8 -*-
#
#  christoffel_calculator.py

#  Copyright 2014 Solkar (Sk)
#  <http://www.quantenforum.de/memberlist.php?mode=viewprofile&u=482>
#
#  version history:
#       0.1 (buggy and widely untested) - Feb 10th 2014 (Sk)   
#            @BUG yet just one possible dump per run
#                 (constraint realized by main() anyway)
#                 likely this a side effect related to sympy symbol allocation
#            @TODO implement simpification of generated symbolics if possible
#            @TODO design a class for encapsulation of the yet globals
#            @TODO add means to attach symbolic indices t,r.. instead of 0,1.. 
#                  to the Γ symbols 
#       0.1.1 (still buggy and widely untested) - Feb 11th 2014 (Sk)   
#            Fixed a typo in FLRW metric, improved Unicode output
#            BUGs and TODOs as before
#       0.1.2 de-Sitter metric included by Bernhard
#
###############################################################################

#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 of the License, or
#  (at your option) any later version.

#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.

#  You should have received a copy of the GNU General Public License
#  along with this program; if not, write to the Free Software
#  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
#  MA 02110-1301, USA.

###############################################################################
"""
    sympy¹-based snippets for generating symbolic representations of
    Christoffel symbols

    GR-related conventions:
    c := 1, t leading, index numbering starts at 0, +2 signature

    ¹http://sympy.org/en/index.html

    Bibliography:
    [MTW73] Misner, C.; Thorne, K. & Wheeler, J. Gravitation.
            W. H. Freeman, 1973.
"""

import sys       # for argv, stderr 
import itertools # for index tuple generation
import codecs
UTF8Writer = codecs.getwriter('utf8')
sys.stdout = UTF8Writer(sys.stdout)

import sympy
import sympy.diffgeom as dg
import sympy.printing.pretty as pp

DIM = 4
CHRISTOFFEL_INDEX_TUPLES = itertools.product(range(DIM), repeat = 3)

LATEX_ENCODER = lambda t,C: '\Gamma_{' + str(t[1]) + str(t[2]) + '}^{' + str(t[0]) + '} = ' + sympy.latex(C)
#UTF8_ENCODER = lambda t,C: 'Γ^' + str(t[0]) + '_' + str(t[1]) + str(t[2]) + " = " + repr(C)
UTF8_ENCODER = lambda t,C: pp(sympy.Symbol(u'\u0393' + "^" + str(t[0]) + '_' + str(t[1]) + str(t[2])),wrap_line = False) + unicode(" =\n") + pp(C)

# + " = "

ENCODER_LABELS=['utf-8','latex']
ENCODERS = {s:f for s,f in zip(ENCODER_LABELS, [UTF8_ENCODER, LATEX_ENCODER])}

def dump_christoffels(christoffels, encoder):
    '''
        delegates output of non-vanishing Γ^i_jk to <encoder>
    '''
    for (i,j,k) in CHRISTOFFEL_INDEX_TUPLES:
        C = christoffels[i][j][k]
        if C != 0:
            print(encoder((i,j,k),C))

   


M = dg.Manifold('M', DIM)
SSP_CS = dg.CoordSystem('Schwarzschild_polar', dg.Patch('P', M),
            ['t', 'r', 'theta', 'phi'])

#
# it's not exactly nice polluting the global namespace with lowercase
# symbols which are also CS-dependent, but python does neither provide
# file-local mods like Perl nor a "using" directive like C++11
#
# @TODO Design a singleton class for encapsulation
#
t, r, theta, phi = SSP_CS.coord_functions() 
dt2, dr2, dtheta2, dphi2 = [dg.TensorProduct(dx,dx)
                                    for dx in SSP_CS.base_oneforms()] 
#

def ds2_FLRW_SSP():
    """
        returns symbolic representation
        of FLRW line element (see e.g. [MTW73] eqns(27.24,27.27))
        in Schwarzschild polar coordinates
    """
    k = sympy.symbols('k')
    dsigma2 = ((1/(1-(k*r**2))) * dr2 + r**2 * (dtheta2 + (sympy.sin(theta))**2 * dphi2))
    a = sympy.symbols('a',cls=sympy.Function)
    ds2 = -dt2 + a(t)**2 * dsigma2
    return ds2
   
def ds2_SS_SSP(use_Rs = True):
    '''
        returns symbolic representation
        of Schwarzschild line element
        in Schwarzschild polar coordinates

        use_rs determines whether "2*M" or "R_s" is used for the symbolics
    '''
    ds2 = None
    f = None
    if use_Rs:
        R_s = sympy.symbols('R_s')
        f = lambda r : 1-R_s/r
    else:
        m = sympy.symbols('m')
        f = lambda r : 1 - 2*m/r
    ds2 = -f(r)*dt2 + 1/f(r) * dr2 + r**2 * (dtheta2 + (sympy.sin(theta))**2 * dphi2)
    return ds2
   
def ds2_deSitter_SC():
    """
        returns symbolic representation
        of deSitter line element (see Wikipedia)
        in static coordinates
    """
    ds2 = None
    a = sympy.symbols('a')
    f1 = lambda r : 1-r**2/a

    ds2 = -f1(r)*dt2 + 1/f1(r) * dr2 + r**2 * (dtheta2 + (sympy.sin(theta))**2 * dphi2)
   
    return ds2

DS_LABELS = ('FLRW', 'Schwarzschild', 'deSitter')
DS_MAKERS = {s:f for s,f in zip(DS_LABELS, [ds2_FLRW_SSP,ds2_SS_SSP,ds2_deSitter_SC])}


def usage():
    s = "usage: python2.7 christoffel_calculator.py ("
    s += "|".join(DS_LABELS)
    s += ") ("
    s += "|".join(ENCODER_LABELS)
    s += ")\ne.g.:"
    s += "\n\tpython2.7 christoffel_calculator.py Schwarzschild utf-8"
    return s


def main():
    """
        @TODO check why just single invocation (as enforced here) runs are
              possible;
              invoking dump_christoffels(..) more than once with
              e.g. a different metric would quietly fail on the #>1 run
        @BUG  dto.
    """
    try:
        # this is MEANT to throw on invalid invocation following
        # "EAFP" paradigma
        # "Easier to Ask for Forgiveness than Permission"
        # see http://docs.python.org/2/glossary.html
        metric_name = sys.argv[1]
        encoding = sys.argv[2]
        ds2 = DS_MAKERS[metric_name]()
        encoder = ENCODERS[encoding]           

        dump_christoffels(dg.metric_to_Christoffel_2nd(ds2), encoder)
       
    except (IndexError, KeyError) as e:
        sys.stderr.write(usage() + "\n")

if __name__ == '__main__':
    main()


Zu beachten: $\alpha^2$ habe ich durch a abgekürzt.
Freundliche Grüße, B.

langlebiges Kaon

Beiträge: 444

Registriert: Mi 26. Okt 2011, 15:34

Beitrag Do 28. Mai 2015, 20:59

Re:

Bernhard hat geschrieben:Hallo Solkar,
ich habe das Skript mal etwas erweitert, so dass man mit dem Befehl "python christoffel_calculator.py deSitter"


Hi Bernhard!

SUPER Sache, das!

Jetzt fehlt ja eigentlich "nur" :mrgreen: noch Kerr-Newman zum Glück.
Wg QS - Hast Du schon mal zu Kerr-Newman christoffelt und einige bekanntermassen korrekte Christoffels dazu oder eine externe Ref mit dto zur Hand?

Beste Grüsse,
Solkar
Old MacDonald had a form
$e_i\,\wedge\,e_i = 0$
(N.N. on Physics Forums)

Anti-Neutron

Beiträge: 921

Registriert: So 16. Jan 2011, 18:11

Beitrag Do 28. Mai 2015, 21:42

Solkar hat geschrieben:Wg QS - Hast Du schon mal zu Kerr-Newman christoffelt und einige bekanntermassen korrekte Christoffels dazu oder eine externe Ref mit dto zur Hand?

Hallo Solkar,

ich habe vor längerer Zeit mal ein selbstgeschriebenes CAS-Tool mit Hilfe der Berechnung der Komponenten des Ricci-Tensors der Kerr-Metrik getestet. Das Programm war da schon einige Zeit beschäftigt und es lief auch erst dann durch, als ich die kovarianten Komponenten der Metrik gemäß diesem Wikipedia-Abschnitt explizit vorgegeben habe. Diese Komponenten stammen aus dem MTW. Kerr-Newman wäre dann die "volle Ausbaustufe" für Schwarze Löcher, für die mir aktuell aber die Motivation eher fehlt.
br
Freundliche Grüße, B.

langlebiges Kaon

Beiträge: 444

Registriert: Mi 26. Okt 2011, 15:34

Beitrag Fr 29. Mai 2015, 16:40

Re: Kleiner SymPy-basierter Christoffel-Rechner für Python

Ich hab mal Kerr-Newman q'n'd einpflegt; wie zu erwarten sind die Terme wahre Bandwürmer.
Ohne QS i.S.v Abgleich gg gesichertermaßen korrekte Terme nuetzt das natürlich noch wenig.

Aber es ist ein Anfang:

  Code:

0.1.3 BEWARE! This is even more buggy than before - May 29th 2015(Sk) 
        started algorithmic decomposition ("init" function) 
        started adding Kerr/Kerr-Newman metrices


  Code:
#!/usr/bin/python2.7
# -*- coding: utf-8 -*-
#
#  christoffel_calculator.py
#
#  Copyright 2014 Solkar (Sk)
#  <http://www.quantenforum.de/memberlist.php?mode=viewprofile&u=482>
#
#  version history:
#       0.1 (buggy and widely untested) - Feb 10th 2014 (Sk)   
#            @BUG yet just one possible dump per run
#                 (constraint realized by main() anyway)
#                 likely this a side effect related to sympy symbol allocation
#            @TODO implement simpification of generated symbolics if possible
#            @TODO design a class for encapsulation of the yet globals
#            @TODO add means to attach symbolic indices t,r.. instead of 0,1..
#                  to the Γ symbols
#       0.1.1 (still buggy and widely untested) - Feb 11th 2014 (Sk)   
#            Fixed a typo in FLRW metric, improved Unicode output
#            BUGs and TODOs as before
#       0.1.2 de-Sitter metric included by Bernhard
#       0.1.3 BEWARE! This is even more buggy than before - May 29th 2015(Sk) 
#             started algorithmic decomposition ("init" function) 
#             started adding Kerr/Kerr-Newman metrices
#             
###############################################################################
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  You should have received a copy of the GNU General Public License
#  along with this program; if not, write to the Free Software
#  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
#  MA 02110-1301, USA.
#
###############################################################################
"""
    sympy¹-based snippets for generating symbolic representations of
    Christoffel symbols

    GR-related conventions:
    G := 1, c := 1, t leading, index numbering starts at 0, +2 signature

    ¹http://sympy.org/en/index.html

    Bibliography:
    [MTW73] Misner, C.; Thorne, K. & Wheeler, J. Gravitation.
            W. H. Freeman, 1973.
    [Wal84] Wald, R. General Relativity. University of Chicago Press, 1984.           
"""

import sys       # for argv, stderr
import itertools # for index tuple generation
import codecs
UTF8Writer = codecs.getwriter('utf8')
sys.stdout = UTF8Writer(sys.stdout)

import sympy
import sympy.diffgeom as dg
import sympy.printing.pretty as pp

DIM = 4
CHRISTOFFEL_INDEX_TUPLES = itertools.product(range(DIM), repeat = 3)

LATEX_ENCODER = lambda t,C: '\Gamma_{' + str(t[1]) + str(t[2]) + '}^{' + str(t[0]) + '} = ' + sympy.latex(C)
#UTF8_ENCODER = lambda t,C: 'Γ^' + str(t[0]) + '_' + str(t[1]) + str(t[2]) + " = " + repr(C)
UTF8_ENCODER = lambda t,C: pp(sympy.Symbol(u'\u0393' + "^" + str(t[0]) + '_' + str(t[1]) + str(t[2])),wrap_line = False) + unicode(" =\n") + pp(C)

# + " = "

ENCODER_LABELS=['utf-8','latex']
ENCODERS = {s:f for s,f in zip(ENCODER_LABELS, [UTF8_ENCODER, LATEX_ENCODER])}

def dump_christoffels(christoffels, encoder):
    '''
        delegates output of non-vanishing Γ^i_jk to <encoder>
    '''
    for (i,j,k) in CHRISTOFFEL_INDEX_TUPLES:
        C = christoffels[i][j][k]
        if C != 0:
            print(encoder((i,j,k),C))

   


M = dg.Manifold('M', DIM)
SSP_CS = dg.CoordSystem('Schwarzschild_polar', dg.Patch('P', M),
            ['t', 'r', 'theta', 'phi'])

BL_CS = dg.CoordSystem('Boyer-Lindquist', dg.Patch('P', M),
            ['t', 'r', 'theta', 'phi'])

#
# it's not exactly nice polluting the global namespace with lowercase
# symbols which are also CS-dependent, but python does neither provide
# file-local mods like Perl nor a "using" directive like C++11
#
# @TODO Design a singleton class for encapsulation
#
t, r, theta, phi = (None for _ in range(4))
dt2, dr2, dtheta2, dphi2, dt, dp, dtp, dpt = (None for _ in range(8))

# redundant, but required to patch around sympy issues

def init(CS):
    global t, r, theta, phi, dt2, dr2, dtheta2, dphi2, dp, dt, dtp, dpt
    t, r, theta, phi = CS.coord_functions()
    dt2, dr2, dtheta2, dphi2 = [dg.TensorProduct(dx,dx)
                                        for dx in CS.base_oneforms()]
    dt = CS.base_oneforms()[0]
    dp = CS.base_oneforms()[3]
    dtp = dg.TensorProduct(dt,dp)

    # redundant, but required to patch around sympy issues
    dpt = dg.TensorProduct(dp,dt)
    #


def ds2_FLRW_SSP():
    """
        returns symbolic representation
        of FLRW line element (see e.g. [MTW73] eqns(27.24,27.27))
        in Schwarzschild polar coordinates
    """
    k = sympy.symbols('k')
    dsigma2 = ((1/(1-(k*r**2))) * dr2 + r**2 * (dtheta2 + (sympy.sin(theta))**2 * dphi2))
    a = sympy.symbols('a',cls=sympy.Function)
    ds2 = -dt2 + a(t)**2 * dsigma2
    return ds2

   
def ds2_SS_SSP(use_Rs = True):
    '''
        returns symbolic representation
        of Schwarzschild line element
        in Schwarzschild polar coordinates

        use_rs determines whether "2*M" or "R_s" is used for the symbolics
    '''
    ds2 = None
    f = None
    if use_Rs:
        R_s = sympy.symbols('R_s')
        f = lambda r : 1-R_s/r
    else:
        m = sympy.symbols('m')
        f = lambda r : 1 - 2*m/r
    ds2 = -f(r)*dt2 + 1/f(r) * dr2 + r**2 * (dtheta2 + (sympy.sin(theta))**2 * dphi2)
    return ds2

   
def ds2_deSitter_SC():
    """
        returns symbolic representation
        of deSitter line element (see Wikipedia)
        in static coordinates
    """
    ds2 = None
    a = sympy.symbols('a')
    f1 = lambda r : 1-r**2/a

    ds2 = -f1(r)*dt2 + 1/f1(r) * dr2 + r**2 * (dtheta2 + (sympy.sin(theta))**2 * dphi2)
   
    return ds2


def ds2_KerrNewman_BL(e=sympy.symbols('e')):
    """
        returns symbolic representation
        of Kerr-Newman line element as of [Wal84] eq (12.3.1)
        in Boyer-Lindquist coordinates
    """
    ds2 = None
    #J = sympy.symbols('J')
    #M = sympy.symbols('M')
    R_s = sympy.symbols('R_s') # = 2M
    a = sympy.symbols('a') #J/M
    Sigma = lambda r, theta: r**2 + a**2 * (sympy.cos(theta))**2
    Delta = lambda r: r**2 + a**2 + e**2 - R_s*r
   
    g_tt = lambda r, theta: -(Delta(r) - a**2 * sympy.sin(theta)**2)/Sigma(r,theta)
    g_tp = lambda r, theta: -(2 * a * sympy.sin(theta)**2 * (r**2 + a**2 - Delta(r)))/Sigma(r,theta)
    g_pp = lambda r, theta: (((r**2 + a**2)**2 - Delta(r) * a**2 * sympy.sin(theta)**2)/Sigma(r,theta)) *  sympy.sin(theta)**2
    g_rr = lambda r, theta: Sigma(r,theta)/Delta(r)
    g_thth = lambda r, theta: Sigma(r,theta)

    ds2 = (g_tt(r,theta)*dt2
           + .5 * g_tp(r,theta)*dtp + .5 * g_tp(r,theta)*dpt # patching around for strange sympy issues
           + g_pp(r,theta)*dphi2 + g_rr(r,theta)*dr2 + g_thth(r,theta)*dtheta2
    )
    return ds2


def ds2_Kerr_BL():
    """
       like ds2_KerrNewman_BL without charge
    """
    return ds2_KerrNewman_BL(e=0)


DS_LABELS = ('FLRW', 'Schwarzschild', 'deSitter', 'Kerr-Newman', 'Kerr')
DS_MAKERS = {s:f for s,f in
    zip(DS_LABELS, [ds2_FLRW_SSP, ds2_SS_SSP, ds2_deSitter_SC, ds2_KerrNewman_BL, ds2_Kerr_BL])
}


def usage():
    s = "usage: python2.7 christoffel_calculator.py ("
    s += "|".join(DS_LABELS)
    s += ") ("
    s += "|".join(ENCODER_LABELS)
    s += ")\ne.g.:"
    s += "\n\tpython2.7 christoffel_calculator.py Schwarzschild utf-8"
    return s


def main():
    """
        @TODO check why just single invocation (as enforced here) runs are
                possible;
                invoking dump_christoffels(..) more than once with
                e.g. a different metric would quietly fail on the #>1 run
        @BUG  dto.
    """
    try:
        # this is MEANT to throw on invalid invocation following
        # "EAFP" paradigma
        # "Easier to Ask for Forgiveness than Permission"
        # see http://docs.python.org/2/glossary.html
        metric_name = sys.argv[1]
        encoding = sys.argv[2]
        init(BL_CS if metric_name in ['Kerr', 'Kerr-Newman'] else SSP_CS)
       
        ds2 = DS_MAKERS[metric_name]()
        encoder = ENCODERS[encoding]           

        dump_christoffels(dg.metric_to_Christoffel_2nd(ds2), encoder)
       
    except (IndexError, KeyError) as e:
        sys.stderr.write(usage() + "\n")

if __name__ == '__main__':
    main()




EDIT 21:00 Solkar
QUARK! So wie ich das versucht hab, kann ich um sympys issues natürlich nicht patchen. So stimmr zwar $ds$ betragsmäßig, aber die Metrik, und somit die Christoffel nicht mehr.
Old MacDonald had a form
$e_i\,\wedge\,e_i = 0$
(N.N. on Physics Forums)

langlebiges Kaon

Beiträge: 444

Registriert: Mi 26. Okt 2011, 15:34

Beitrag Fr 29. Mai 2015, 19:56

Re: Kleiner SymPy-basierter Christoffel-Rechner für Python

In Version 0.1.3 muss ich also den Code für Kerr/KerrNewman noch bugfixen; das taugt so gar nix.
Old MacDonald had a form
$e_i\,\wedge\,e_i = 0$
(N.N. on Physics Forums)

langlebiges Kaon

Beiträge: 444

Registriert: Mi 26. Okt 2011, 15:34

Beitrag Sa 30. Mai 2015, 16:30

Re: Kleiner SymPy-basierter Christoffel-Rechner für Python

Solkar hat geschrieben:In Version 0.1.3 muss ich also den Code für Kerr/KerrNewman noch bugfixen; das taugt so gar nix.

EDIT 30. May 2015 17:30 Das mag zwar sein, aber nicht aus dem Grund, den ich meinte.
Old MacDonald had a form
$e_i\,\wedge\,e_i = 0$
(N.N. on Physics Forums)

langlebiges Kaon

Beiträge: 444

Registriert: Mi 26. Okt 2011, 15:34

Beitrag Fr 26. Feb 2016, 03:39

Re: Kleiner SymPy-basierter Christoffel-Rechner für Python

Neue Version 0.2.0


  Code:
#!/usr/bin/python3
# -*- coding: utf-8 -*-
#
#  christoffel_calculator.py
#
#  Copyright 2014 Solkar (Sk)
#  <http://www.quantenforum.de/memberlist.php?mode=viewprofile&u=482>
#
#  version history:
#       0.1 (buggy and widely untested) - Feb 10th 2014 (Sk)
#            @BUG yet just one possible dump per run
#                 (constraint realized by main() anyway)
#                 likely this a side effect related to sympy symbol allocation
#            @TODO implement simpification of generated symbolics if possible
#            @TODO design a class for encapsulation of the yet globals
#            @TODO add means to attach symbolic indices t,r.. instead of 0,1..
#                  to the Γ symbols
#       0.1.1 (still buggy and widely untested) - Feb 11th 2014 (Sk)
#            Fixed a typo in FLRW metric, improved Unicode output
#            BUGs and TODOs as before
#       0.1.2 de-Sitter metric included by Bernhard
#       0.1.3 BEWARE! This is even more buggy than before - May 29th 2015(Sk)
#             started algorithmic decomposition ("init" function)
#             started adding Kerr/Kerr-Newman metrices
#                        -----
#       0.2.0 Major overhaul - Feb 26th 2016 (Sk)
#             - python3 ready (python3.4)
#             - PEP8 compliant (skip E731 for lambda - expressions)
#             - added basic expression simplification
#             - added support for a Levi-Civita cylinder metric
#
###############################################################################
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  You should have received a copy of the GNU General Public License
#  along with this program; if not, write to the Free Software
#  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
#  MA 02110-1301, USA.
#
###############################################################################
"""
    sympy¹-based snippets for generating symbolic representations of
    Christoffel symbols

    GR-related conventions:
    G := 1, c := 1, t leading, index numbering starts at 0, +2 signature

    ¹http://sympy.org/en/index.html

    Bibliography:
    [GHW10] Gralla, S. E.; Harte, A. I. & Wald, R. M.
            Bobbing and kicks in electromagnetism and gravity.
            prd, 2010, 81, 104012.
            doi:10.1103/PhysRevD.81.104012
            http://arxiv.org/abs/1004.0679
    [MTW73] Misner, C.; Thorne, K. & Wheeler, J. Gravitation.
            W. H. Freeman, 1973.
    [Wal84] Wald, R. General Relativity. University of Chicago Press, 1984.
"""

import sys        # for argv, stderr
import itertools  # for index tuple generation

import sympy
from sympy import cancel, powsimp, trigsimp
import sympy.diffgeom as dg
import sympy.printing.pretty as pp

DIM = 4
CHRISTOFFEL_INDEX_TUPLES = itertools.product(range(DIM), repeat=3)


class LaTeXEncoder(object):
    def __call__(self, t, C):
        return pp(
            '\Gamma_{' + str(t[1]) + str(t[2]) + '}^{' +
            str(t[0]) + '} = ' + sympy.latex(C)
        )


class UTF8Encoder(object):
    def __call__(self, t, C):
        return pp(
            sympy.Symbol(u'\u0393' + "^" + str(t[0]) + '_' +
                         str(t[1]) + str(t[2])), wrap_line=False
        ) + u' =\n' + pp(C)


ENCODERS = {
    'utf-8': UTF8Encoder(),
    'latex': LaTeXEncoder()
}


def dump_christoffels(christoffels, encoder):
    '''
        delegates output of non-vanishing Γ^i_jk to <encoder>
    '''
    for (i, j, k) in CHRISTOFFEL_INDEX_TUPLES:
        C = christoffels[i][j][k]
        if C != 0:
            # @TODO: Improve simplification
            for f in (trigsimp, powsimp):
                C = f(C)
            print(encoder((i, j, k), C))


M = dg.Manifold('M', DIM)

SSP_CS = dg.CoordSystem('Schwarzschild_polar', dg.Patch('P', M),
                        ['t', 'r', 'theta', 'phi'])

BL_CS = dg.CoordSystem('Boyer-Lindquist', dg.Patch('P', M),
                       ['t', 'r', 'theta', 'phi'])

CYL_CS = dg.CoordSystem('Cylindrical', dg.Patch('P', M),
                        ['t', 'r', 'phi', 'z'])


#
# it's not exactly nice polluting the global namespace with lowercase
# symbols which are also CS-dependent, but python does neither provide
# file-local mods like Perl nor a "using" directive like C++11
#
# @TODO Design a singleton class for encapsulation
#
t, r, theta, phi, z = (None for _ in range(5))
dt2, dr2, dtheta2, dphi2, dz2, dt, dp, dtp, dpt = (None for _ in range(9))


def init(CS):
    global t, r, theta, phi, z
    global dt2, dr2, dtheta2, dphi2, dz2, dp, dt, dtp, dpt

    if CS == CYL_CS:
        t, r, phi, z = CS.coord_functions()
        dt2, dr2, dphi2, dz2 = [dg.TensorProduct(dx, dx)
                                for dx in CS.base_oneforms()]
    else:
        t, r, theta, phi = CS.coord_functions()
        dt2, dr2, dtheta2, dphi2 = [dg.TensorProduct(dx, dx)
                                    for dx in CS.base_oneforms()]

    dt = CS.base_oneforms()[0]
    dp = CS.base_oneforms()[3]
    dtp = dg.TensorProduct(dt, dp)

    # redundant, but required to patch around sympy issues
    dpt = dg.TensorProduct(dp, dt)
    #


def ds2_FLRW_SSP():
    """
        returns symbolic representation
        of FLRW line element (see e.g. [MTW73] eqns(27.24,27.27))
        in Schwarzschild polar coordinates
    """
    k = sympy.symbols('k')
    dsigma2 = ((1/(1-(k*r**2))) * dr2 +
               r**2 * (dtheta2 + (sympy.sin(theta))**2 * dphi2))
    a = sympy.symbols('a', cls=sympy.Function)
    ds2 = -dt2 + a(t)**2 * dsigma2
    return ds2


def ds2_SS_SSP(use_Rs=True):
    '''
        returns symbolic representation
        of Schwarzschild line element
        in Schwarzschild polar coordinates

        use_rs determines whether "2*M" or "R_s" is used for the symbolics
    '''
    f = None
    if use_Rs:
        R_s = sympy.symbols('R_s')
        f = lambda r: 1-R_s/r
    else:
        m = sympy.symbols('m')
        f = lambda r: 1 - 2*m/r

    ds2 = (
        -f(r)*dt2 + 1/f(r) * dr2 +
        r**2 * (dtheta2 + (sympy.sin(theta))**2 * dphi2)
    )
    return ds2


def ds2_deSitter_SC():
    """
        returns symbolic representation
        of deSitter line element (see Wikipedia)
        in static coordinates
    """
    ds2 = None
    a = sympy.symbols('a')
    f1 = lambda r: 1-r**2/a

    ds2 = (
        -f1(r)*dt2 + 1/f1(r) * dr2 +
        r**2 * (dtheta2 + (sympy.sin(theta))**2 * dphi2)
    )

    return ds2


def ds2_KerrNewman_BL(e=sympy.symbols('e')):
    """
        returns symbolic representation
        of Kerr-Newman line element as of [Wal84] eq (12.3.1)
        in Boyer-Lindquist coordinates
    """
    R_s = sympy.symbols('R_s')  # = 2M
    a = sympy.symbols('a')  # J/M
    Sigma = lambda r, theta: r**2 + a**2 * (sympy.cos(theta))**2
    Delta = lambda r: r**2 + a**2 + e**2 - R_s*r

    g_tt = lambda r, theta: (
        -(Delta(r) - a**2 * sympy.sin(theta)**2)/Sigma(r, theta)
    )
    g_tp = lambda r, theta: (
        -(2 * a * sympy.sin(theta)**2 * (r**2 + a**2 - Delta(r))) *
        1/Sigma(r, theta)
    )
    g_pp = lambda r, theta: (
        (
            ((r**2 + a**2)**2 - Delta(r) *
             a**2 * sympy.sin(theta)**2)/Sigma(r, theta)
        ) * sympy.sin(theta)**2
    )
    g_rr = lambda r, theta: Sigma(r, theta)/Delta(r)
    g_thth = lambda r, theta: Sigma(r, theta)

    ds2 = (g_tt(r, theta)*dt2 +
           # patching around for strange sympy issues
           .5 * g_tp(r, theta)*dtp + .5 * g_tp(r, theta)*dpt +
           #
           g_pp(r, theta)*dphi2 + g_rr(r, theta)*dr2 +
           g_thth(r, theta)*dtheta2)

    return ds2


def ds2_Kerr_BL():
    """
       like ds2_KerrNewman_BL without charge
    """
    return ds2_KerrNewman_BL(e=0)


def ds2_LeviCivita_CYL():
    """
        returns symbolic representation
        of line element for [GHW10] eq (91)
        in cylindrical coordinates
    """
    ds2 = None
    rho = sympy.symbols('rho')
    alpha = sympy.symbols('alpha')

    ds2 = (- r**(4*rho) * dt2 +
           r**(4*rho(2*rho - 1)) * (dr2 + dz2) +
           alpha**(-2) * r**(2-4*rho) * dphi2)

    return ds2


DS_MAKERS = {
    'FLRW':          ds2_FLRW_SSP,
    'Schwarzschild': ds2_SS_SSP,
    'deSitter':      ds2_deSitter_SC,
    'Kerr-Newman':   ds2_KerrNewman_BL,
    'Kerr':          ds2_Kerr_BL,
    'Levi-Civita':   ds2_LeviCivita_CYL
}


def usage():
    s = "usage: python3 christoffel_calculator.py ("
    s += "|".join(DS_MAKERS.values())
    s += ") ("
    s += "|".join(ENCODERS.keys())
    s += ")\ne.g.:"
    s += "\n\tpython3 christoffel_calculator.py Schwarzschild utf-8"
    return s


def main():
    """
        @TODO check why just single invocation (as enforced here) runs are
                possible;
                invoking dump_christoffels(..) more than once with
                e.g. a different metric would quietly fail on the #>1 run
        @BUG  dto.
    """
    try:
        # this is MEANT to throw on invalid invocation following
        # "EAFP" paradigma
        # "Easier to Ask for Forgiveness than Permission"
        # see http://docs.python.org/2/glossary.html
        metric_name = sys.argv[1]
        encoding = sys.argv[2]
        CS = {
            'Kerr': BL_CS,
            'Kerr-Newman': BL_CS,
            'FLRW': SSP_CS,
            'Schwarzschild': SSP_CS,
            'deSitter': SSP_CS,
            'Levi-Civita': CYL_CS
        }[metric_name]

        init(CS)

        ds2 = DS_MAKERS[metric_name]()
        encoder = ENCODERS[encoding]

        dump_christoffels(dg.metric_to_Christoffel_2nd(ds2), encoder)

    except (IndexError, KeyError) as e:
        sys.stderr.write(usage() + "\n")

if __name__ == '__main__':
    main()
Old MacDonald had a form
$e_i\,\wedge\,e_i = 0$
(N.N. on Physics Forums)

Zurück zu Relativitätstheorie und Kosmologie

Wer ist online?

Mitglieder in diesem Forum: 0 Mitglieder und 5 Gäste

cron
Powered by phpBB® Forum Software © phpBB Group
Designed by ST Software for PTF.
Deutsche Übersetzung durch phpBB.de