EMANの物理学 過去ログ No.6832 〜

 ● テンソル/リーマン幾何学の計算ソフト

  投稿者:小林@那須 - 2009/05/01(Fri) 10:01  No.6832 
小林@那須と申します。python sf と呼んでいる計算ソフトを作っています。python の scipy, sympy, vpython, pylab と幾つかの自作モジュールを組み合わせてプリプロセッサを介して計算させます。下にあるような計算をさせられます。傲慢と取られるのは分かっていますが、大部分のユーザーに対してMatlab/Mathematica を超えられたと自負しています。

http://www.nasuinfo.or.jp/FreeSpace/kenji/sf/sfShotConcise/sfShotConcise.htm

汎用的に使える Jacobian を今回実装できたので、自分が一般相対論を理解することも兼ねて、共変微分の数値計算に使ってみました。付録にあるコードのように Christoffel 接続テンソルを含んだ計算を簡単に実行できます。

このテンソル計算について皆様の意見を伺わせてください。

○ python sf ではテンソルの足と縮約インデックスを明記することを、できるだけ避けています。デフォルトのテンソル積では、テンソルの右端と左端のインデックスどうしを縮約させています。行列積の拡張になっています。付録コードの共変微分の計算では、これで済ませられました。

テンソルの縮約インデックスを明記して操作すると、どうしても複雑になってしまいます。C や Java で書いたのと殆ど同じになってしまいます。この意味で縮約インデックスを省略できることは、テンソル計算させるソフトとして必須の要求仕様だと思っています。

でも縮約でのインデックスを明記するテンソル積操作に慣れた方々に受け入れられるのか心配です。如何でしょうか?

○ コンピュータによるテンソル計算の状況を google で検索してみたのですが、インデックスを明示的に使った、複雑な形式のものしか見つけられませんでした。でも探し方がまずいのかもしれません。既にテンソル計算を簡潔に行わせている方、そのような例を知っている方がいらしたら教えてもらえますでしょうか。

=========== 付録コード begin ==================
∂J(...) で Jacobian を計算させています。
//@@
# http://homepage2.nifty.com/eman/relativity/co_dif.html
# 座標変換後の共変微分値が微分値の座標変換値と一致することを数値確認する
def fx(X,Y):
  return ~[X^2+0.1Y, 2Y^2]

def fX(x,y):
  return ~[(x-0.1(0.5y)**0.5)**0.5, (0.5*y)**0.5 ]

G__ = ~[[1,0],
    [0,1]]

pos = (1.5, 2)
p_xy = fx(*pos)
print "x:",fx(*pos)
print "X:",fX(*p_xy)

fJxy_JXY = λ X,Y: ~[
  # ∂x`[i]
  # ----------
  #   ∂X`[j]
  [ 2X, 0.1],
  [ 0, 4Y]
  #[ 2X, 0.0],
  #[ 0, 4Y]
          ]

fJXY_Jxy = λ x,y: ~[
  # inverted and transposed matrix function
  # ∂X`[i]
  # ----------
  #   ∂x`[j]
  [0.5*(x - 0.0707106781186548*y**0.5)**-0.5,
  -0.0176776695296637*y**-0.5*(x - 0.0707106781186548*y**0.5)**-0.5],
  #[0.5*x**-0.5, 0],
  [ 0, 0.353553390593274*y**-0.5]
          ]

# X,Y 座標系で A = ~[3,4] のベクトルが空間に一様に分布しているとする
# この A を x,y 座標系で表すと、下のように位置に依存して変化する量になる
def A_(X,Y):
  #return ~[3, 4]
  #return ~[3X^2, 4Y^2+1]
  return ~[3X^2+Y, 4Y^2+1]
def a_xy(x,y):
  #     ∂X[j]
  #  A[j] ----------
  #      ∂x`[i]
  return A_(*fX(x,y)) fJXY_Jxy(x,y)

def fg__(x,y):
  mtAt = fJXY_Jxy(x,y)
  return mtAt.t G__ mtAt

print "********************"
print "Γ`__ ≡ fJxy_JXY(*pos) ∂J(fJXY_Jxy)(*p_xy)"
Γ`__ = fJxy_JXY(*pos) ∂J(fJXY_Jxy)(*p_xy)
print Γ`__
print " ∇ (a_xy) - a_xy Γ`__ ≡ Jacobian(a_xy) - a_xy Γ`__"
print ∂J(a_xy)(*p_xy) - a_xy(*p_xy) Γ`__

print "********************"
print "ベクトル A_ の∇:Jacobian 値"
print ∂J(A_)(*pos)
# 下の .t:transpose に気がつかなくて惑わされた
print "Jacobian(A_) を x,y 座標に変換した行列値が上の共変微分値に一致する"
print fJXY_Jxy(*p_xy).t ∂J(A_)(*pos) fJXY_Jxy(*p_xy)
//@@@
//copy \#####.### temp.py /y
//testPy.py -fs temp.py
//__tempConverted.py
=========== 付録コード end ==================
x: [ 2.45 8. ]
---- ClTensor ----
X: [ 1.5 2. ]
---- ClTensor ----
********************
Γ`__ ≡ fJxy_JXY(*pos) ∂J(fJXY_Jxy)(*p_xy)
[[[ -2.22222222e-01 2.77777778e-03]
[ 2.77777778e-03 -3.47222222e-05]]

[[ 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 -6.25000000e-02]]]
---- ClTensor ----
∇ (a_xy) - a_xy Γ`__ ≡ Jacobian(a_xy) - a_xy Γ`__
[[ 1. 0.02916667]
[-0.0125 0.24963542]]
---- ClTensor ----
********************
ベクトル A_ の∇:Jacobian 値
[[ 9. 1.]
[ 0. 16.]]
---- ClTensor ----
Jacobian(A_) を x,y 座標に変換した行列値が上の共変微分値に一致する
[[ 1. 0.02916667]
[-0.0125 0.24963542]]
---- ClTensor ----

  投稿者:小林@那須 - 2009/05/19(Tue) 07:10  No.6898  <Home>
python sf を使って、EMAN さんの Web page 第5部「リーマン幾何学」をトレースしてみました。なんとかリーマン幾何・一般相対論がイメージできるようになってきました。現在、一般相対論を次のように理解しています。

・アインシュタイン・テンソル  ${G^{i,j}}$  と運動量テンソル  ${T^{i,j}}$  が正比例する。
・ ${G~{i,j}}$  はリーマン幾何学の Ricci tensor と  ${R^{i,j}}$  とよく似ている。 ${-0.5 g^{i,j} R}$  が違うだけ。
ただし R はRicci scalar
・Ricci tensor は、二次元、三次元では Rieman tensor と実質的に同じもの。
・Rieman tensor は空間の回転変形をあらわすもの。
<== Rieman tensor は、Christoffel 接続係数と その空間微分の組み合わせで表される
<== Rieman tensor は、計量 g`` の二階までの空間微分の組み合わせでも表される

短く言うと、運動量テンソル  ${T^{i,j}}$  が、アインシュタイン方程式によって計量  ${g^{i,j}}$  の分布を 1 から変化させる。即ち空間を歪める。ただし、この方程式から空間の歪み:計量  ${g^{i,j}}$  を計算するのは非常に難しい。一般解は見つかっていない。

一般相対性理論を理解しようとするとき、Christoffel 接続係数や計量から Rieman 曲率/Ricci テンソルを導く過程に多くの手間が係ります。三次元を超える多次元空間の話であり、幾何学でありながら、幾何学的なイメージを構築できるようになるのは普通の人間には不可能に近い話に思えます。その過程は膨大なテンソル式の羅列を理解していく作業になるのですが、それは三つ/四つもあるテンソルの足を使っての縮約の繰り返しです。各テンソル式の幾何学的意味を理解しながらのテンソル式を変形していく作業は無理難題に思えます。途中で誤りが簡単に入り込みます。リーマン幾何学から一般相対論に至る道程の最初の段階で、何をやっているのか分からなくなってしまうのが普通の人間だと思います。

EMAN さんの web page の説明は、幾何学的イメージの構築に拘って説明を続けています。自然・数学法則を理解することを目的とした説明で一貫しています。抽象数学による一般化・厳密化された説明とは対極にあるものです。数式の写経だけでは納得できない EMAN さんの真骨頂が出ていると思います。普通の人間に、リーマン幾何学の数式を使った幾何学的イメージを掴ませるのに最も成功している web ページだと思います。この web ページがなかったら、私の能力では最初の Christoffel 接続係数の段階で幾何学的イメージを構築することができずに、脱落していたと思います。

今は  ${\Gamma^{i}_{j,k}}$ :Christoffel 接続係数を次のように理解しています。誤りがあったら指摘してやってください。

「Christoffel 接続係数は本質的には  ${\partial{^2}X/\partial{x} \partial{x}}$  です。x-->X 座標変換のベクトル値関数を二回 Jacobian 微分したものです。これは空間の歪みと見なせます。単純な一変数の scalar 関数の場合だと  ${\dif{^2f}{x^2}}$ :二階微分に相当します。ただし  ${\dif{^2f}{x^2}}$  は scalar 値ですが、 ${\partial{^2}X/\partial{x} \partial{x}}$  はテンソルです。しかも歪んだ空間でのテンソルです。反変、共変どちらで見るのかを指定する必要があります。 この意味で、Christoffel 接続係数  ${Gamma^{i}_{j,k}}$  は  ${\partial{^2}X/\partial{x} \partial{x}}$  三階テンソルの前に、 ${\pdif{x}{X}}$  二階テンソルが掛かっています。歪んでいない Euclid 空間 X からの変換行列が掛かっています。これにより、歪んだ空間での共変ベクトル v の空間微分を求めるとき、歪んだ空間の寄与分を  ${v_{i } \Gamma^{i}_{ j,k}}$  即ち v * Γ`__ と表せます。」


  投稿者:小林@那須 - 2009/05/19(Tue) 07:14  No.6899  <Home>
EMAN さんの web page が幾何学イメージを説明していると先に書きました。しかし私の能力では、 EMAN さんの説明だけでは まだ幾何学イメージを構成していけませんでした。Python sf で、Christoffel テンソルが計量テンソルと、それらの Jacobian を組み合わせた数値計算を、電卓でのように行わせていくことで、やっと幾何学的イメージが掴めてきました。Python sf を使ったコンピュータの助けがなかったら、Christoffel 接続係数の段階で振り落とされたと思います。

Christoffel/計量テンソルを電卓での様に計算させるため、Euclid 空間:X から歪んだ空間 x への座標変換ベクトル関数を [X^2+0.1Y, 2Y^2] と仮定し、これを元に下のようなリーマン幾何学で使う変数・関数を sfCrrnt.py ファイルに用意してやりました。これにより宣言することなく、下の関数変数を使えるようになります。

-------------- sfCrrnt.py で定義してある変数・関数 begin ----------------
fx, fX 二次元座標変換ベクトル関数
   X_ が Euclid 座標から、歪んだ空間 x_ への変換関数
   具体的には、下のベクトル値を返す関数です。x_ は下の関数です。
    def fx(X,Y):
      return ~[X^2+0.1Y, 2Y^2]

      J(fx)
fJxX(X,Y) ------: Jacobian 微分だが、解析解の関数:
      J(X)

      J(fX)
fJXx(x,y) ------: Jacobian 微分だが、解析解の関数
      J(x)

    fJxX, fJXx は数値解ではなく、解析解を使って書いてあります。
    ∵計算誤差が入り込みやすい微分計算を減らすためです。

pX   Euxclid 空間における位置。[2,3] が入れてある
px   歪んだ空間における pX に対応する位置。[ 1.2, 8]が入っている

Γ`__  Christoffel coeffcient at [2,3]

fgx__  g_{i,j} matrix function
fgx   g^{i,j} matrix function

A_(X,Y) Euclidian 座標系でのベクトル分布。~[3* X**2+Y, 4* Y**2+1] のベクトル分布としています。

a_xy(x,y): A_ の分布を歪んだ空転での分布として表したもの。
     <== A_(*X_(x,y)) fJXx(x,y) を返しているだけ

a_ Euclidian 座標系の [2,3] 位置におけるベクトル A_ の x_ 座標によるベクトル値
g__ Euclidian 座標系の [2,3] 位置における、座標変換の微分から定まる計量
g`` Euclidian 座標系の [2,3] 位置における計量逆行列

jXx Euclidian 座標系の [2,3] 位置における X -> x 座標変換の Jacobian 値
jxX Euclidian 座標系の [2,3] 位置における x -> X 座標変換の Jacobian 値

∇v_:共変ベクトルの共変微分関数(a_xy(x,y) などのベクトル分布関数を引数にとる)
∇m__:二階の共変テンソル:即ち行列の共変微分関数(fgx__(x,y) など行列分布関数を引数にとる)

∂JX ∂J と同じ Jacobian 碑文だが、∂/∂X を明示するため導入する
∂Jx ∂J と同じ Jacobian 碑文だが、∂/∂x_ を明示するため導入する
-------------- RmnEx で定義してある変数・関数 end ----------------

これらの関数・変数を前提にできると、多くの計算が以下のようにワンライナーで可能になります。テンソル積の計算で、インデックスを明記せずに行列のときと同様に計算させられるときはワンライナーで書けることが多くなります。

g__ g``
===============================
[[ 1.00000000e+00 2.22044605e-15]
[ 4.44089210e-16 1.00000000e+00]]
---- ClTensor ----


g__ g`` の計算では数値微分は入っておらず、解析解の組み合わせ計算だけですが、対角項は 0 ではなく、10^-16 のオーダーの誤差が入ってきます。位置を Jacobian 微分したときの関数 fJXx は 16 桁の浮動少数点で書かれているからです。

計量 g__ の共変微分は、その分布関数 fgx__ に二階の共変テンソルの微分関数 ∇m__ を適用して計算ます。
∇m__(fgx__)(*px)
===============================
[[[ -6.35206193e-11 5.94767291e-14]
[ -3.80894795e-12 -1.23883651e-14]]

[[ -3.80894773e-12 -1.23883651e-14]
[ -2.25513076e-13 -1.06559293e-13]]]
---- ClTensor ----

理論どおり 0 になります。∇m__(...) は関数を返しているので、px 以外の場所、たとえば x == (1.5,1.8) での g__ の共変微分値も同様に計算できます。他の場所でも計量 g__ の共変微分は 0 になります。

∇m__(fgx__)(1.5, 1.8)
===============================
[[[ -1.84270738e-10 -1.43999938e-14]
[ -5.64926578e-13 1.22020585e-14]]

[[ -5.64926198e-13 1.22020449e-14]
[ 3.97515273e-14 -4.49890195e-11]]]
---- ClTensor ----

数値微分なので、どうしても誤差が付きまといます。でも Richardson 補外法を使っているので、この場合は四桁の悪化で済んでいます。二回の Jacobian 微分までなら、多くの場合六桁の精度まで期待できるようです。数式を誤って理解していないかの確認には十分に使えます。

なお、Euclid 空間から歪んだ空間への変換が [X^2+0.1Y, 2Y^2] の関数ベクトルのときは、Rieman tensor が 0 になってしまいます。一般的に [f(x), g(x)] のように R^2 --> R^2 への一対一の変換関数で空間を歪めただけでは Rieman tensor は 0 になってしまうようです。球面のように R^2 --> R^3 のように、異なる次元への変換ベクトル関数を用意してやれば Rieman tensor が 0 でなくなります。でも R^2--> R^2 への変換では、その Jacobian に逆行列が存在します。Christoffel 接続係数の性質を調べるときには、こちらの方が便利です。


以下に、sfCrrnt.py の具体的な中身を示します。

----------------- sfCrrnt.py の中身 begin ------------------
from __future__ import division
# -*- encoding: cp932 -*-
"""' Rieman 幾何学を数値シミュレーションするときの関数群、座標変換関数 x_, X_,
 その Jacobian: fJxX, fJXx, ベクトル分布 A_, a_, 計量 fgXx__ を定義する。
 位置: pX=(1,2), px = x_(*pX)
 Chirstoffel 分布関数: fΓ`__(x,y)
 位置 px におけるChirstoffel 係数値: Γ`__
 共変微分 ∇v_, ∇m__

'"""
from __future__ import division
# -*- encoding: cp932 -*-
from sf import *
#from sfCrrntIni import *
def fx(X,Y):
    #return krry([X**2+0.1* Y, 2* Y**2])
    #return krry([X**2+0.1* X* Y**2, 2*Y**2])
    return krry([X**2+0.1* X* Y**3, 2*Y**2])

def fX(x,y):
    #return krry([(x-0.1* (0.5* y)**0.5)**0.5, (0.5*y)**0.5 ])
    #return krry([-0.025*y + 1/2*(4*x + 0.0025*y**2)**(1/2)
    #, 2**(1/2)*y**(1/2)/2
    return krry([ 1/2*(4*x + 0.00125*(y**(3/2))**2)**(1/2) - 0.0125*2**(1/2)*y**(3/2)
            , 2**(1/2)*y**(1/2)/2
            ])

G__ = krry([[1,0],
        [0,1]])

pX = (1,2)
px = fx(*pX)

fJxX = lambda X,Y: krry(
    # ∂x`[i]
    # ----------
    # ∂X`[j]
    #[[ 2* X, 0.1],
    # [ 0, 4* Y]

    #[[ 2* X+0.1* Y**2, 0.2* X* Y]
    #,[ 0, 4* Y]
    [[ 2* X+0.1* Y**3, 0.3* X* Y**2]
    ,[ 0, 4* Y]
    ])

fJXx = lambda x,y: krry(
    # inverted and transpXed matrix function
    # ∂X`[i]
    # ----------
    # ∂x`[j]
    #[[0.5*(x - 0.0707106781186548*y**0.5)**-0.5,
    # -0.0176776695296637*y**-0.5*(x - 0.0707106781186548*y**0.5)**-0.5],
    # [ 0, 0.353553390593274*y**-0.5]
    [[(4*x + 0.00125*(y**1.5)**2)**-0.5
                  , -0.0265165042944955*y**0.5 + 0.0009375*y**-1.*(y**1.5)**2*(4*x
                   + 0.00125*(y**1.5)**2)**-0.5]
    ,[ 0, 0.353553390593274*y**-0.5]
    ])
                    
def fgx__(x,y):
    return fJXx(x,y).t * G__ * fJXx(x,y)

def k_fgx_bq__bq____(x,y):
    return fJxX(*fX(x,y)) * fJxX(*fX(x,y)).t

g__=fgx__(*px)
k_g_bq__bq____=k_fgx_bq__bq____(*px)

jXx = fJXx(*px)
jxX = fJxX(*pX)

def A_(X,Y):
    return krry([3* X**2+Y, 4* Y**2+1])

def a_xy(x,y):
    # ∂X[j]
    # A[j] ----------
    # ∂x`[i]
    return A_(*fX(x,y)) * fJXx(x,y)

a_ = a_xy(*px)

k_f_lGamma__bq______ = lambda x,y:fJxX(*fX(x,y)) * Jc_(fJXx)(x,y)
k__lGamma__bq______ = k_f_lGamma__bq______(*px)

def k__Nabla_v____(v_):
    return lambda x,y:Jc_(v_)(x,y) - v_(x,y) * fJxX(*fX(x,y)) * Jc_(fJXx)(x,y)

def k__Nabla_m_____(m__):
    def innerF(x,y):
        k__lGamma__bq______ = k_f_lGamma__bq______(x,y)
        tns = Jc_(m__)(x,y) - m__(x,y) * k__lGamma__bq______
        for i,j,k in mrng(2,2,2):
            tns[i,j,k] -= k__lGamma__bq______[:,i,k] * m__(x,y)[:,j]

        return tns

    return innerF

k__Round_Jx___ = Jc_
k__Round_JX___ = Jc_

"""'
//@@
def fx(X,Y):
    return ~[X^2+0.1Y, 2Y^2]

def fX(x,y):
    return ~[(x-0.1(0.5y)**0.5)**0.5, (0.5*y)**0.5 ]

G__ = ~[[1,0],
        [0,1]]

pX = (1,2)
px = fx(*pX)

fJxX = λ X,Y: ~[
    # ∂x`[i]
    # ----------
    # ∂X`[j]
    [ 2X, 0.1],
    [ 0, 4Y]
                     ]

fJXx = λ x,y: ~[
    # inverted and transpXed matrix function
    # ∂X`[i]
    # ----------
    # ∂x`[j]
    [0.5*(x - 0.0707106781186548*y**0.5)**-0.5,
      -0.0176776695296637*y**-0.5*(x - 0.0707106781186548*y**0.5)**-0.5],
    [ 0, 0.353553390593274*y**-0.5]
                     ]
                    
def fgx__(x,y):
    return fJXx(x,y).t G__ fJXx(x,y)

def fgx``(x,y):
    return fJxX(*fX(x,y)) fJxX(*fX(x,y)).t

g__=fgx__(*px)
g``=fgx``(*px)

jXx = fJXx(*px)
jxX = fJxX(*pX)

def A_(X,Y):
    return ~[3X^2+Y, 4Y^2+1]

def a_xy(x,y):
    # ∂X[j]
    # A[j] ----------
    # ∂x`[i]
    return A_(*fX(x,y)) fJXx(x,y)

a_ = a_xy(*px)

fΓ`__ = λ x,y:fJxX(*fX(x,y)) Jc_(fJXx)(x,y)
Γ`__ = fΓ`__(*px)

def ∇v_(v_):
    return λ x,y:Jc_(v_)(x,y) - v_(x,y) fJxX(*fX(x,y)) Jc_(fJXx)(x,y)

def ∇m__(m__):
    def innerF(x,y):
        Γ`__ = fΓ`__(x,y)
        tns = Jc_(m__)(x,y) - m__(x,y) Γ`__
        for i,j,k in mrng(2,2,2):
            tns[i,j,k] -= Γ`__[:,i,k] m__(x,y)[:,j]

        return tns

    return innerF

∂Jx = Jc_
∂JX = Jc_
//@@@
//copy \#####.### temp.py /y
//testPy.py -fs temp.py
//type __tempConverted.py


'"""
----------------- sfCrrnt.py の中身 end ------------------

  投稿者:小林@那須 - 2009/05/19(Tue) 07:27  No.6900  <Home>
上のような sfCrrnt.py を前提に Christoffel 接続係数と計量の関係を調べていきます。

下側インデックスだけの Γ___:Christoffel 係数を座標変換関数の Jacobian および Jacobian Jacobian の積として計算します。通常の足を書き連ねたテンソル計算よりも、こちらの方が幾何学的な意味が分かりやすいと思うのですが、いかがてしょうか。足の記号の交換よりも行列の transpose のほうが違いが明確に表現されます。

Γ___ at pX;; fJXx(*px).t ∂Jx(fJXx)(*px)
===============================
[[[-0.03253853 -0.00195231]
 [-0.00195231 0.00012202]]

[[ 0.00488078 0.00029285]
 [ 0.00029285 -0.00099487]]]
---- ClTensor ----

g`` Γ___ を計算;; g`` fJXx(*px).t ∂Jx(fJXx)(*px)
===============================
[[[ -2.55102041e-01 -1.53061224e-02]
 [ -1.53061225e-02 -8.41836735e-03]]

[[ -5.25980297e-16 -3.15588178e-17]
 [ -3.15588178e-17 -6.25000000e-02]]]
---- ClTensor ----

これは Γ`__ すなわち fJxX(*X_(x,y)) ∂J(fJXx)(x,y) と同じになります。(ちなみに球面では この関係式が成り立たなくなります。)
Γ`__;; fJxX(*pX) ∂Jx(fJXx)(*px)
===============================
[[[-0.25510204 -0.01530612]
 [-0.01530612 -0.00841837]]

[[ 0. 0. ]
 [ 0. -0.0625 ]]]
---- ClTensor ----

共変微分の章:http://homepage2.nifty.com/eman/relativity/co_dif.html の (1) 式は、次のように計算できます。
左辺値 ∂g__[i,j]/∂x[k];;∂Jx(fgx__)(*px)
===============================
[[[-0.06507705 -0.00390462]
 [ 0.00292847 0.00041487]]

[[ 0.00292847 0.00041487]
 [ 0.00058569 -0.00198973]]]
---- ClTensor ----

(1) 式の右辺値は、足が入り組んでおり、右端と左端の足の縮約にできないので、ワンライナーで書くのは無理があります。下のようにループを回して計算します。上の g__ の Jacobian 微分に一致します。

//@@
JJXxAt = ∂Jx(fJXx)(*px)

tns = kzrs(2,2,2)
for i,j,k in mrng(2,2,2):
tns[i,j,k] = JJXxAt[:,i,k] JXx[:,j] + JXx[:,i] JJXxAt[:,j,k]

pp(tns)
print nearlyEq(tns, ∂Jx(fgx__)(*px)) # g__ の Jacobian と 10^-6 の精度で同じことを確認
//@@@
//copy \#####.### temp.py /y
//testPy.py -fs temp.py
//__tempConverted.py

[[[ -0.0650771, -0.00390462]
,[ 0.00292847, 0.000414866]]
,[[ 0.00292847, 0.000414866]
,[ 0.000585693, -0.00198973]]]
True

このループ計算は、C, Java, Matlab などでは四重のループ構文で書かねばなりません。python sf では一重のループで済ませられ、数式の記述に近いことが分かってもらえるかと思います。可読性も十分にあります。

計量 fgx__ の変化率:Jacobian によって Γ___, Γ`__ を計算させます。∂^X/∂x ∂x から計算したものと同じ値になります。ここらの計算をワンライナーで行わせるか、複数行で行わせるかの判断は、人によって分かれると思います。

jgx=∂Jx(fgx__)(*px);tns=kzrs(2,2,2);for k,i,j in mrng(2,2,2):tns[k,i,j]=0.5(jgx[j,k,i]+jgx[k,i,j] - jgx[i,j,k]);pp(tns)
[[[ -0.0325385, -0.00195231]
,[ -0.00195231, 0.000122019]]
,[[ 0.00488078, 0.000292847]
,[ 0.000292847, -0.000994865]]]

jgx=∂Jx(fgx__)(*px);tns=kzrs(2,2,2);for k,i,j in mrng(2,2,2):tns[k,i,j]=0.5 g``[k,:] (jgx[j,:,i]+jgx[:,i,j] - jgx[i,j,:]);pp(tns)
[[[ -0.255102, -0.0153061]
,[ -0.0153061, -0.00841837]]
,[[ 0, 0]
,[ 0, -0.0625]]]


Christoffel 係数が空間の接平面の凸凹具合 ∂^2 X/∂x ∂x に座標変換を掛けたものであること、それが計量の位置変化すなわち Jacobian:∂J(fgx):(これも空間の凸凹具合) と一次変換で関係付けるられることが分かります。

通常、Christoffel 接続係数まわりの計算を C/Java なでのプログラムで行おうとすると、プログラムのデバッグ作業で大きく時間をとられてしまいます。物理をやりたいはずなのに、プログラム・デバッグ作業ばかりをしていることになりがちです。Matlab や Mathematica でも、それに近いものがあります。

特定の事例のときに限られるとはいえ、Christoffel 接続係数に関連す計算を、上のようにワンライナーでできてしまうので物理に集中できます。プログラムを組むというより、電卓操作に近いからです。物理や数学の式の意味を数値計算で確認するときに行いたいのは、ここでやった電卓でのような計算です。何十行にも渡るプログラムの作成による計算ではありません。この意味で多くのユーザーにとって、python sf は Matlab や Mathematica を超えたと主張します。如何でしょうか?


----------------------------------------------------
「平行移動」の章には式変形の途中に誤りがあります。これは、誤り報告として別スレッドにします。

  投稿者:小林@那須 - 2009/05/25(Mon) 17:00  No.6924  <Home>
Ricci tensor を計量  ${g_{ij}}$  と その微分で表してみました。

sympy では  ${\partial{^2}f/\partial{x}\partial{y} - \partial{^2}f/\partial{y}\partial{x}}$ が 0 になってませんでした。下の URL にある、このバグの対策パッチを作ってくれたので、Ricci tensor を計量  ${g_{ij}}$  と その微分で表してみました。

http://code.google.com/p/sympy/issues/detail?id=1435

「場の古典論」(92.4) 式に次のような公式がありました。

<tex>R_{imkl} = 1/2(\frac{\partial{^2}g_{im}}{\partial{x_{k}}\partial{x_{l}}}              +\frac{\partial{^2}g_{kl}}{\partial{x_{i}}\partial{x_{m}}}              -\frac{\partial{^2}g_{im}}{\partial{x_{k}}\partial{x_{m}}}              -\frac{\partial{^2}g_{km}}{\partial{x_{i}}\partial{x_{l}}}) \\                        + g_{np}(\Gamma^{n}_{ kl} \Gamma^{p}_{ im} - \Gamma^{n}_{ km} \Gamma^{p}_{ il})</tex>

なので ${g_{np}(\Gamma^{n}_{ kl} \Gamma^{p}_{ im} - \Gamma^{n}_{ km} \Gamma^{p}_{ il})}$ の項だけを計算しました。

下のようになりました。後ろ側の項を  ${g_{ij}}$  で表現すると、予想以上に ごちゃごちゃになります。あたりまえですが、 ${g_{ij}}$  が分母に入り込んできます。Ricci tensor を計量 ${g_{ij}}$ で表した数式が教科書に 出てこないのが納得できました。

下のコードは python と 先のパッチをあてた sympy064 が入っていれば動くようにしてあります。興味のある方は試してみてください。

----------- code begin -----------------
//@@
import sympy as ts

N=2

# multiple range generator
def mrng(*args):
    """' multiple range generator by Kenji Kobayashi in kVerifeir Lab
        for sf:scientific functional calculator
                http://www.nasuinfo.or.jp/FreeSpace/kenji/index.htm

        Kenji has the copy rights of this code

        list(mrng(2,3))"
        ===============================
        [(0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (1, 2)]

        list(mrng((1,4),(5,7)))
        ===============================
        [(1, 5), (1, 6), (2, 5), (2, 6), (3, 5), (3, 6)]

        list(mrng((1,6,2),(5,9,3)))

        ===============================
        [(1, 5), (1, 8), (3, 5), (3, 8), (5, 5), (5, 8)]

        an example of usage
            dctAt={}
            for i,j in mrng(10,20):
                dctAt[i,j] = i+j
               
    '"""
    head, tail = args[0], args[1:]
    if type(head) in [int, long, float]:
        head = range(head)
    elif isinstance(head, tuple) or isinstance(head, list):
        head = range(*head)
    if tail:
        for i in head:
            for j in mrng(*tail):
                if isinstance(j, tuple):
                    yield (i,)+j
                else:
                    yield (i, j)
    else:
        for i in head:
            yield i

x,y = ts.symbols("xy")
lstSmbl=[x,y]

g__ = ts.Matrix( N,N, lambda i,j:ts.Function('g__'+str(i)+str(j))(x,y) )
guu = g__.inv()
#print guu

Gamma_ddd = {}
for i,j,k in mrng(N,N,N):
    Gamma_ddd[i,j,k] = ts.Rational(1,2)*(g__[i,j].diff(lstSmbl[k])
                                       + g__[k,i].diff(lstSmbl[j])
                                       - g__[j,k].diff(lstSmbl[i])
                                        )
Gamma_udd = {}
for i,j,k in mrng(N,N,N):
    Gamma_udd[i,j,k] = sum([guu[i,m]* Gamma_ddd[m,j,k] for m in range(N)])

R____tail = {}
for i,k,l,m in mrng(N,N,N,N):
    R____tail[i,k,l,m] = sum([
            g__[n,p]*( Gamma_udd[n,k,l]*Gamma_udd[p,i,m]
                     - Gamma_udd[n,k,m]*Gamma_udd[p,i,l] )
                                for n,p in mrng(N,N)
                             ])

#print R____tail

R__tail = {}
for i,j in mrng(N,N):
    R__tail[i,j] = sum([ R____tail[m,i,m,j] for m in range(N)])

#print R__tail[0, 0].expand()
#print "======================"
print R__tail[0, 0]
//@@@
----------- code end -----------------

result

(((D(g__01(x, y), y)/2 + D(g__11(x, y), x)/2 - D(g__10(x, y), y)/2)/(-g__01(x,
y)*g__10(x, y)/g__00(x, y) + g__11(x, y)) - (D(g__00(x, y), y)/2 + D(g__01(x,
y), x)/2 - D(g__10(x, y), x)/2)*g__10(x, y)/((-g__01(x, y)*g__10(x,
y)/g__00(x, y) + g__11(x, y))*g__00(x, y)))*((D(g__10(x, y), y)/2 + D(g__11(x,
y), x)/2 - D(g__01(x, y), y)/2)/(-g__01(x, y)*g__10(x, y)/g__00(x, y) +
g__11(x, y)) - (D(g__00(x, y), y)/2 + D(g__10(x, y), x)/2 - D(g__01(x, y),
x)/2)*g__10(x, y)/((-g__01(x, y)*g__10(x, y)/g__00(x, y) + g__11(x,
y))*g__00(x, y))) - ((D(g__01(x, y), x)/2 + D(g__10(x, y), x)/2 - D(g__00(x,
y), y)/2)/(-g__01(x, y)*g__10(x, y)/g__00(x, y) + g__11(x, y)) - D(g__00(x,
y), x)*g__10(x, y)/(2*(-g__01(x, y)*g__10(x, y)/g__00(x, y) + g__11(x,
y))*g__00(x, y)))*(D(g__11(x, y), y)/(2*(-g__01(x, y)*g__10(x, y)/g__00(x, y)
+ g__11(x, y))) - (D(g__01(x, y), y)/2 + D(g__10(x, y), y)/2 - D(g__11(x, y),
x)/2)*g__10(x, y)/((-g__01(x, y)*g__10(x, y)/g__00(x, y) + g__11(x,
y))*g__00(x, y))))*g__11(x, y) + (((D(g__10(x, y), y)/2 + D(g__11(x, y), x)/2 -
D(g__01(x, y), y)/2)/(-g__01(x, y)*g__10(x, y)/g__00(x, y) + g__11(x, y)) -
(D(g__00(x, y), y)/2 + D(g__10(x, y), x)/2 - D(g__01(x, y), x)/2)*g__10(x,
y)/((-g__01(x, y)*g__10(x, y)/g__00(x, y) + g__11(x, y))*g__00(x,
y)))*((g__01(x, y)*g__10(x, y)/((-g__01(x, y)*g__10(x, y)/g__00(x, y) +
g__11(x, y))*g__00(x, y)**2) + 1/g__00(x, y))*(D(g__00(x, y), y)/2 +
D(g__01(x, y), x)/2 - D(g__10(x, y), x)/2) - (D(g__01(x, y), y)/2 + D(g__11(x,
y), x)/2 - D(g__10(x, y), y)/2)*g__01(x, y)/((-g__01(x, y)*g__10(x,
y)/g__00(x, y) + g__11(x, y))*g__00(x, y))) - ((D(g__01(x, y), x)/2 +
D(g__10(x, y), x)/2 - D(g__00(x, y), y)/2)/(-g__01(x, y)*g__10(x, y)/g__00(x,
y) + g__11(x, y)) - D(g__00(x, y), x)*g__10(x, y)/(2*(-g__01(x, y)*g__10(x,
y)/g__00(x, y) + g__11(x, y))*g__00(x, y)))*((g__01(x, y)*g__10(x,
y)/((-g__01(x, y)*g__10(x, y)/g__00(x, y) + g__11(x, y))*g__00(x, y)**2) +
1/g__00(x, y))*(D(g__01(x, y), y)/2 + D(g__10(x, y), y)/2 - D(g__11(x, y),
x)/2) - D(g__11(x, y), y)*g__01(x, y)/(2*(-g__01(x, y)*g__10(x, y)/g__00(x, y)
+ g__11(x, y))*g__00(x, y))))*g__10(x, y) + (((D(g__01(x, y), y)/2 +
D(g__11(x, y), x)/2 - D(g__10(x, y), y)/2)/(-g__01(x, y)*g__10(x, y)/g__00(x,
y) + g__11(x, y)) - (D(g__00(x, y), y)/2 + D(g__01(x, y), x)/2 - D(g__10(x,
y), x)/2)*g__10(x, y)/((-g__01(x, y)*g__10(x, y)/g__00(x, y) + g__11(x,
y))*g__00(x, y)))*((g__01(x, y)*g__10(x, y)/((-g__01(x, y)*g__10(x,
y)/g__00(x, y) + g__11(x, y))*g__00(x, y)**2) + 1/g__00(x, y))*(D(g__00(x, y),
y)/2 + D(g__10(x, y), x)/2 - D(g__01(x, y), x)/2) - (D(g__10(x, y), y)/2 +
D(g__11(x, y), x)/2 - D(g__01(x, y), y)/2)*g__01(x, y)/((-g__01(x, y)*g__10(x,
y)/g__00(x, y) + g__11(x, y))*g__00(x, y))) - ((g__01(x, y)*g__10(x,
y)/((-g__01(x, y)*g__10(x, y)/g__00(x, y) + g__11(x, y))*g__00(x, y)**2) +
1/g__00(x, y))*D(g__00(x, y), x)/2 - (D(g__01(x, y), x)/2 + D(g__10(x, y),
x)/2 - D(g__00(x, y), y)/2)*g__01(x, y)/((-g__01(x, y)*g__10(x, y)/g__00(x, y)
+ g__11(x, y))*g__00(x, y)))*(D(g__11(x, y), y)/(2*(-g__01(x, y)*g__10(x,
y)/g__00(x, y) + g__11(x, y))) - (D(g__01(x, y), y)/2 + D(g__10(x, y), y)/2 -
D(g__11(x, y), x)/2)*g__10(x, y)/((-g__01(x, y)*g__10(x, y)/g__00(x, y) +
g__11(x, y))*g__00(x, y))))*g__01(x, y) + (((g__01(x, y)*g__10(x,
y)/((-g__01(x, y)*g__10(x, y)/g__00(x, y) + g__11(x, y))*g__00(x, y)**2) +
1/g__00(x, y))*(D(g__00(x, y), y)/2 + D(g__10(x, y), x)/2 - D(g__01(x, y),
x)/2) - (D(g__10(x, y), y)/2 + D(g__11(x, y), x)/2 - D(g__01(x, y),
y)/2)*g__01(x, y)/((-g__01(x, y)*g__10(x, y)/g__00(x, y) + g__11(x,
y))*g__00(x, y)))*((g__01(x, y)*g__10(x, y)/((-g__01(x, y)*g__10(x,
y)/g__00(x, y) + g__11(x, y))*g__00(x, y)**2) + 1/g__00(x, y))*(D(g__00(x, y),
y)/2 + D(g__01(x, y), x)/2 - D(g__10(x, y), x)/2) - (D(g__01(x, y), y)/2 +
D(g__11(x, y), x)/2 - D(g__10(x, y), y)/2)*g__01(x, y)/((-g__01(x, y)*g__10(x,
y)/g__00(x, y) + g__11(x, y))*g__00(x, y))) - ((g__01(x, y)*g__10(x,
y)/((-g__01(x, y)*g__10(x, y)/g__00(x, y) + g__11(x, y))*g__00(x, y)**2) +
1/g__00(x, y))*D(g__00(x, y), x)/2 - (D(g__01(x, y), x)/2 + D(g__10(x, y),
x)/2 - D(g__00(x, y), y)/2)*g__01(x, y)/((-g__01(x, y)*g__10(x, y)/g__00(x, y)
+ g__11(x, y))*g__00(x, y)))*((g__01(x, y)*g__10(x, y)/((-g__01(x, y)*g__10(x,
y)/g__00(x, y) + g__11(x, y))*g__00(x, y)**2) + 1/g__00(x, y))*(D(g__01(x, y),
y)/2 + D(g__10(x, y), y)/2 - D(g__11(x, y), x)/2) - D(g__11(x, y), y)*g__01(x,
y)/(2*(-g__01(x, y)*g__10(x, y)/g__00(x, y) + g__11(x, y))*g__00(x,
y))))*g__00(x, y)