見出し画像

一般相対論、曲率テンソルの計算プログラム(Python)

物理になりたい人。

面倒な曲率テンソルを計算するためのpythonプログラム、書きました。
何を計算しているかが分かりやすいようにプログラムを書いたつもりです。（添え字に関する対称性を使えばもう少し早くなるかもしれないが、計算オーダーは変わらないため、わかりやすさを優先しました。）

ファイルのダウンロードは↓から。


曲率等の定義

・Christoffel 記号：

$$
\Gamma^\mu\ _{\rho\sigma} = \frac{1}{2} g^{\mu \nu}(\partial_{\sigma} g_{\nu\rho} + \partial_{\rho} g_{\sigma\nu} - \partial_{\nu} g_{\rho\sigma})
$$

・Riemannテンソル

$$
R^\alpha\, _{\beta\mu\nu} = \partial_\mu \Gamma^\alpha\, {\nu\beta} - \partial_\nu \Gamma^\alpha\,_{\mu\beta} + \Gamma^\alpha\,_{\mu\tau}\Gamma^\tau\,_{\nu\beta} - \Gamma^\alpha\,_{\nu\tau} \Gamma^\tau\,_{\mu\beta}
$$

・Ricciテンソル

$$
R_{\mu \nu} = R^\alpha\,_{\mu\alpha\nu}
$$

・Ricciスカラー

$$
R = g^{\mu \nu} R_{\mu \nu}
$$

・Einsteinテンソル

$$
G_{\mu \nu} = R_{\mu \nu} - \frac{1}{2} g_{\mu \nu} R
$$

プログラム

関数等の定義

import sympy

class curvature:
    """
    計量テンソルから各種曲率を計算する.
    曲率の定義：
    ・Christoffel symbols :
        \Gamma^\mu\,_{\rho\sigma} = \frac{1}{2} g^{\mu \nu}(\partial_{\sigma} g_{\nu\rho} + \partial_{\rho} g_{\sigma\nu} - \partial_{\nu} g_{\rho\sigma})
    ・Riemann tensor :
        R^\alpha\, _{\beta\mu\nu} = \partial_\mu \Gamma^\alpha\, {\nu\beta} - \partial_\nu \Gamma^\alpha\,_{\mu\beta} + \Gamma^\alpha\,_{\mu\tau}\Gamma^\tau\,_{\nu\beta} - \Gamma^\alpha\,_{\nu\tau} \Gamma^\tau\,_{\mu\beta}
    ・Ricci tensor :
        R_{\mu \nu} = R^\alpha\,_{\mu\alpha\nu}
    ・Scalar curvature :
        R = g^{\mu \nu} R_{\mu \nu}
    ・Einstein tensor :
        G_{\mu \nu} = R_{\mu \nu} - \frac{1}{2} g_{\mu \nu} R
    """
    def __init__(self, vars, g):
        """
        vars : 時空座標変数のリスト(次元は問わない)
        g    : 計量テンソル
        """
        self.xmu = sympy.Matrix([[v] for v in vars])
        self.dim = len(self.xmu)
        if g.shape !=(self.dim, self.dim):
            raise ValueError("計量テンソルの次元と時空次元が一致しません.")
        self.g = g
        self.ginv = g.inv()
        print("ds**2 = "+str(self.ds()))

        # データ初期化
        self.Gamma = None
        self.Riemann = None
        self.Ricci = None
        self.Ricci_scalar = None
        self.Einstein = None

    def ds(self):
        """
        ds**2を式として返す.
        """
        dxmu = sympy.symbols(["d"+str(v) for v in self.xmu])
        ds = 0
        for mu in range(self.dim):
            for nu in range(self.dim):
                ds = ds + dxmu[mu] * dxmu[nu]*self.g[mu,nu]
        return ds

    def cal_Gamma(self):
        """
        Christoffel記号を計算する.
        Γ^a_b_cが self.Gamma[a,b,c] に格納される.
        """
        self.Gamma = sympy.MutableDenseNDimArray.zeros(c.dim,c.dim,c.dim)
        for mu in range(c.dim):
            for rho in range(c.dim):
                for sigma in range(c.dim):
                    gamma = 0
                    for nu in range(self.dim):
                            gamma = gamma + sympy.Rational(1,2)*self.ginv[mu,nu]*(self.g[nu,rho].diff(self.xmu[sigma])+self.g[sigma,nu].diff(self.xmu[rho])-self.g[rho,sigma].diff(self.xmu[nu]))
                    self.Gamma[mu,rho,sigma] = sympy.simplify(gamma)

    def cal_Riemann(self):
        """
        Riemannテンソルを計算する.
        R^a_b_c_dが self.Riemann[a,b,c,d] に格納される.
        """
        if self.Gamma is None:
            self.cal_Gamma()
        self.Riemann = sympy.MutableDenseNDimArray.zeros(c.dim,c.dim,c.dim,c.dim)
        for alpha in range(c.dim):
            for beta in range(c.dim):
                for mu in range(c.dim):
                    for nu in range(c.dim):
                        riemann = self.Gamma[alpha,nu,beta].diff(self.xmu[mu]) - self.Gamma[alpha,mu,beta].diff(self.xmu[nu])
                        for tau in range(self.dim):
                            riemann = riemann + self.Gamma[alpha,mu,tau]*self.Gamma[tau,nu,beta] - self.Gamma[alpha,nu,tau]*self.Gamma[tau,mu,beta]
                        self.Riemann[alpha,beta,mu,nu] = sympy.simplify(riemann)

    def cal_Ricci(self):
        """
        Ricciテンソルを計算する.
        R_a_bが self.Ricci[a,b] に格納される.
        """
        if self.Riemann is None:
            self.cal_Riemann()
        self.Ricci = sympy.MutableDenseNDimArray.zeros(c.dim,c.dim)
        for mu in range(c.dim):
            for nu in range(c.dim):
                ricci = 0
                for alpha in range(self.dim):
                    ricci = ricci + self.Riemann[alpha, mu, alpha, nu]
                self.Ricci[mu,nu] = sympy.simplify(ricci)

    def cal_Ricci_scalar(self):
        """
        Ricciスカラーを計算する.
        計算結果は self.Ricci_scalar に格納される.
        """
        if self.Ricci is None:
            self.cal_Ricci()
        R = 0
        for mu in range(self.dim):
            for nu in range(self.dim):
                R = R + self.ginv[mu,nu]*self.Ricci[mu,nu]
        self.Ricci_scalar = sympy.simplify(R)

    def cal_Einstein(self):
        """
        Einsteinテンソルを計算する.
        G_a_bが self.Einstein[a,b] に格納される.
        """
        if self.Ricci_scalar is None:
            self.cal_Ricci_scalar()
        self.Einstein = sympy.MutableDenseNDimArray.zeros(c.dim,c.dim)
        for mu in range(c.dim):
            for nu in range(c.dim):
                G = self.Ricci[mu,nu] - sympy.Rational(1,2)*self.g[mu,nu]*self.Ricci_scalar
                self.Einstein[mu,nu] = sympy.simplify(G)

使用例

ts_var = ["t","r","theta","phi"]    # 時空変数を文字で指定
sympy.var(ts_var)                   # 上で指定された文字をpython変数に変換
f = sympy.Function("f")(t, r)       # 任意函数を定義

g = sympy.Matrix([
[f, 0, 0, 0],
[0,-1/f, 0, 0],
[0, 0,-r**2, 0],
[0, 0, 0, -r**2*sympy.sin(theta)**2]])# 計量テンソルを定義.

c = curvature(ts_var,g)
c.cal_Gamma()               # Christoffel記号を計算.
print(c.Gamma[1,1,0])       # Gamma^1_1_0を出力.
c.cal_Riemann()             # Riemannテンソルを計算.
print(c.Riemann[1,2,2,1])   # R^1_2_2_1を出力.
c.cal_Ricci()               # Ricciテンソルを計算.
print(c.Ricci[0,1])         # R_0_1を出力.
c.cal_Ricci_scalar()        # Ricciスカラーを計算.
print(c.Ricci_scalar)       # Ricciスカラーを出力.
c.cal_Einstein()            # Einsteinテンソルを計算.
print(c.Einstein[1,1])      # G_1_1を出力.

おまけ：計算結果を整形してtexに出力するプログラム

def output_as_tex(c,output_list = ["g","Christoffel","Riemann","Ricci","Ricci scalar","Einstein"],nonzeroonly=True,filename="out"):
    """
    output_list⊂["g","Christoffel","Riemann","Ricci","Ricci scalar","Einstein"] : 出力する内容を選択.
    nonzeroonly   : Trueを指定すると、各計算結果の非ゼロ成分だけ出力される.
    filename      : 出力texファイルの名前を指定.
    """
    b = "\\begin{align}\n"
    e = "\n\\end{align}\n"
    outtext = r"""
The definitions of curvatures:
\begin{itemize}
    \item Christoffel symbols
        \begin{align}
            \Gamma^\mu\,_{\rho\sigma}
            = \frac{1}{2} g^{\mu \nu}(\partial_{\sigma} g_{\nu\rho}
            + \partial_{\rho} g_{\sigma\nu}
            - \partial_{\nu} g_{\rho\sigma})
        \end{align}

    \item Riemann tensor
        \begin{align}
            R^\alpha\,_{\beta\mu\nu}
            = \partial_\mu \Gamma^\alpha\, {\nu\beta}
            - \partial_\nu \Gamma^\alpha\,_{\mu\beta}
            + \Gamma^\alpha\,_{\mu\tau}\Gamma^\tau\,_{\nu\beta}
            - \Gamma^\alpha\,_{\nu\tau} \Gamma^\tau\,_{\mu\beta}
        \end{align}

    \item Ricci tensor
        \begin{align}
            R_{\mu \nu} = R^\alpha\,_{\mu\alpha\nu}
        \end{align}

    \item Scalar curvature
        \begin{align}
            R = g^{\mu \nu} R_{\mu \nu}
        \end{align}

    \item Einstein tensor
        \begin{align}
        G_{\mu \nu} = R_{\mu \nu} - \frac{1}{2} g_{\mu \nu} R
        \end{align}
\end{itemize}

\hrulefill \\
Calculation condition:
\begin{itemize}
"""
    outtext = outtext + "\t\item spacetime dimension: " + str(c.dim) + "\\\\"
    outtext = outtext + "\t\item spacetime variables: " + str([v for v in c.xmu]) + "\\\\"
    outtext = outtext + "\t\item metric :"
    outtext = outtext + b + "\tds^2 = " + sympy.latex(c.ds()) + e
    outtext = outtext + r"""
\end{itemize}
\hrulefill \\
"""
    if nonzeroonly:
        outtext = outtext + "Calculation result (Only non-zero components)\n"
    else:
        outtext = outtext + "Calculation result\n"

    temptext =""
    if "g" in output_list:
        temptext = temptext + "\\item metric tensor \n"
        temptext = temptext + b + "\t\tg_{\\mu\\nu} = " + sympy.latex(c.g) + e
    if "Christoffel" in output_list:
        temptext = temptext + "\\item Christoffel symbols \n"
        ccount = 0
        for mu in range(c.dim):
            for rho in range(c.dim):
                for sigma in range(c.dim):
                    g = c.Gamma[mu,rho,sigma]
                    if g==0 and nonzeroonly:
                        pass
                    else:
                        temptext = temptext + b + f"\t\t\\Gamma^{mu}\,_{{{rho}{sigma}}} = " + sympy.latex(g) + e
                        ccount +=1
        if nonzeroonly and (ccount == 0):
            temptext = temptext + b + "\t\t\\Gamma^\\mu\,_{\\rho\\sigma} = 0" + e
    if "Riemann" in output_list:
        temptext = temptext + "\\item Riemann tensor \n"
        ccount = 0
        for alpha in range(c.dim):
            for beta in range(c.dim):
                for mu in range(c.dim):
                    for nu in range(c.dim):
                        r = c.Riemann[alpha,beta,mu,nu]
                        if r==0 and nonzeroonly:
                            pass
                        else:
                            temptext = temptext + b + f"\t\tR^{alpha}\,_{{{beta}{mu}{nu}}} = " + sympy.latex(r) + e
                            ccount +=1
        if nonzeroonly and (ccount == 0):
            temptext = temptext + b + "\t\tR^\\alpha\,_{\\beta\\mu\\nu} = 0" + e
    if "Ricci" in output_list:
        temptext = temptext + "\\item Ricci tensor \n"
        ccount = 0
        for mu in range(c.dim):
            for nu in range(c.dim):
                ric = c.Ricci[mu,nu]
                if ric==0 and nonzeroonly:
                    pass
                else:
                    temptext = temptext + b + f"\t\tR_{{{mu}{nu}}} = " + sympy.latex(ric) + e
                    ccount +=1
        if nonzeroonly and (ccount == 0):
            temptext = temptext + b + "\t\tR_\\mu\\nu = 0" + e
    if "Ricci scalar" in output_list:
        temptext = temptext + "\\item Ricci scalar \n"
        temptext = temptext + b + "\t\tR = " + sympy.latex(c.Ricci_scalar) + e
    if "Einstein" in output_list:
        temptext = temptext + "\\item Einstein tensor \n"
        ccount = 0
        for mu in range(c.dim):
            for nu in range(c.dim):
                g = c.Einstein[mu,nu]
                if g==0 and nonzeroonly:
                    pass
                else:
                    temptext = temptext + b + f"\t\tG_{{{mu}{nu}}} = " + sympy.latex(g) + e
                    ccount +=1
        if nonzeroonly and (ccount == 0):
            temptext = temptext + b + "\t\tG_\\mu\\nu = 0" + e
    if temptext != "":
        outtext = outtext + "\\begin{itemize}\n" + temptext + "\n\\end{itemize}"

    with open(f"{filename}.tex","wt") as f:
        f.write(outtext)

↓使用例

output_as_tex(c,nonzeroonly=True,filename="out") # 計算結果をtexファイルに出力


いいなと思ったら応援しよう！