[CAS-Lab] ブラックホールをあなたの手の上に　一般相対性理論２

本稿は下記物理学ノート/一般相対性理論２で数式検証、グラフ描画に使用したMathematicaのプログラムの公開（個人での利用を許可）と簡単な解説となる。

一般相対性理論１ [一般座標変換編]から続いての解説となる。曲がった時空を扱うためにテンソルやリーマン幾何等の数学知識が必要になる。数式を追いかけるだけでは実感がわかないだろうから数式処理 (CAS) を使って視覚化を行い、あなたの理解を助けられるようにしようと思う。

リーマン幾何学で記述された、曲がった空間でのベクトルの動きや接続、曲率などの数式は難解で、その意味するところを直観的に理解することが難しいと思う。本稿に沿って、自分で処理を確認しながら実行することにより、その数式の意味を実感できるようになると思う。

サンプルプログラム実行前の準備

1) Mathematicaの環境構築

Mathematicaの環境が無い場合はこちらを参考に環境を整備してください。

2) サンプルプログラム使用パッケージ（描画関数や座標変換関数等）のセーブ

以下ボタンで［.nb］ファイルをダウンロードしMathematicaで開いて実行してみてください。実行後いくつかのパッケージ（xxx.pac) がWindowsの場合デフォルトで”ドキュメント”　ディレクトリにセーブされるはずです。

(*Mathematica*)

(*---------------基本パッケージ-----------------*)
　　：
　関数、変数定義　省略
　　：
(*---------------変換パッケージ-----------------*)
　　：
　関数、変数定義　省略
　　：
(*---------------描画パッケージ-----------------*)
　　：
　関数、変数定義　省略
　　：
(*---------ベクトル平行移動関連の共通関数--------*)
　　：
　関数、変数定義　省略
　　：
Save[(*cnd*)"udVectTranslate.pac",
 {
 セーブする関数、変数名のリスト  
 }];

上記ではソースコードの転記と解説を省略しました。ダウンロードした実物で内容は確認してください。

サンプルプログラムの実行

以下のボタンでサンプルプログラムをダウンロードしMathematicaで開いて実行できます。

最初にグローバル変数等の初期化と上記でセーブしたパッケージの内、”udPlotPac1.pac” 、 “udElemntPac1.pac” 、”udTransform1.pac”と”udVectTranslate.pac”を読み込ます。これにより以降のプログラムで使用可能となります。

(*初期化、ユーザー定義関数パッケージの読み込み*)
(*Mathematica*)
ClearAll["Global`*"]

Get["udPlotPac1.pac"];
Get["udElemntPac1.pac"];
Get["udTransform1.pac"];
Get["udVectTranslate.pac"];

prON = False;

サンプルプログラムの実行

以上で作業準備は終わりました。以降『一般相対性理論２』の数式検証とグラフ化に関するプログラムになります。各章節事にリンクを貼ってあるので内容を読み比べながらプログラムを実行してください。尚、内容はほぼ同じだが数式や画像の構成に多少の違いがあることをご承知ねがいたい。

2.1 共変微分

反変ベクトルの平行移動

手順１

(*速度ベクトルの移動と共変・反変ベクトルの移動例*)
(*共変ベクトル・反変ベクトルの例*)
prON = False;
(*描画画像の設定*)
P1 = {2, 1}; P2 = {1, 3}; P3 = {1, 1};
opL3 = {P1, P2, P3, P1};
opT3 = {{"b", {0.8, 2}}, {"a", {1.5, 0.8}}, {"c", {1.6, 2.2}}, {"P1", {2.3, 0.7}}, {"P2", {1.2, 3.26}}};
goh = Graphics[{AbsoluteThickness[3], RGBColor[1, 0, 0, 1],
        Line[opL3],
        Table[
     Text[Style[opT3[[ix, 1]], Large, FontFamily -> "Times"], 
      opT3[[ix, 2]]], {ix, Length[opT3]}]
    }];
opT3 = { {"b", {0.8, 2}},  {"a", {1.5, 0.8}}, { 
    "c", {1.6, 2.2}}, {"P1", {2.1, 0.7}}, {"P2", {0.7, 3.0}}};
opTu3 = Table[{opT3[[ix, 1]], X2Up[opT3[[ix, 2]], Be]}, {ix, 
    Length[opT3]}];
(*座標描画領域の設定*)
Gh = {{-1, 3}, {-1, 4}};
Gt = {{0, 4}, {-1, 4}};
(*移動計算用微小係数 表示サイズVdxに掛け、1ループでの移動距離を計算する*)
d = 0.001;
(*地縛参照用反変ベクトルAの表示情報*)
dxA = {0.5, 0.};
(*地縛参照用速度ベクトルVの表示情報*)
dxV = {-0.2, 0.4};
(*共変ベクトル計算用のスカラー場: E *)
　(*EのX系での分布関数*)
SFx = (0.2 (x[[1]] + 0 x[[2]]))^1;
　(*EのU系での分布関数： X\[Rule]U変換で生成*)
SFu = SFx /. Thread[x -> U2Xp[u, Be]];
(*反変・共変ベクトルの選択描画設定関数*)
intW[TorB_] := Module[{Wu, VecColer},
   If[TorB, Wu = (PuAe[[1]] - PuAb[[1]]); VecColer = Blue,
                    Wu = d (PuBe[[1]] - PuAb[[1]]); VecColer = Orange];
   {Wu, VecColer, TorB}];
(*地縛ベクトルの位置*)
PxAb = {{2, 1}, {1.5, 2}, {1, 3}};
PxBb = PxAb;

(*X系視点の表示*)
PxVe = Table[PxAb[[ix]] + d dxV, {ix, Length[PxAb]}];
PxAe = Table[PxAb[[ix]] + d dxA, {ix, Length[PxAb]}];
dxB = {D[SFx, x[[1]]], D[SFx, x[[2]]]}
PxBe = Table[(PxBb[[ix]] + dxB), {ix, Length[PxBb]}];
PuBb = Table[X2Upf[PxBb[[ix]] // N, Be], {ix, Length[PxBb]}];
duB = {D[SFu, u[[1]]], D[SFu, u[[2]]]};
PuBe = Table[(PuBb[[ix]] + duB) /. Thread[u -> PuBb[[ix]]], {ix, Length[PuBb]}];
grVecVx = Graphics[{Thick, RGBColor[0, 0.5, 0.1, 1], Arrow[Transpose[{PxAb, PxAb + (PxVe - PxAb)/d}]]}];
grVecAx = Graphics[{Thick, Blue, Arrow[Transpose[{PxAb, PxAb + (PxAe - PxAb)/d}]]}];
grVecBx = Graphics[{Thick, Orange, Arrow[Transpose[{PxBb, PxBe}]]}];
gxH = Show[grdMesh2D[X2X, Be, {4 {-1, 1}, {-1, 4}}, cBluH, 1]];
gu2xH = Show[grdMesh2D[U2Xp, U2Xp, {4 {0, 1}, 1 Pi {0, 1}}, cGrnH, 100]];
Show[gxH, gu2xH, goh, grVecVx, grVecAx, grVecBx, 
  PlotRange -> {{Gh[[1, 1]] + 0.7, Gh[[1, 2]] + 0.01}, {Gh[[2, 1]] + 0.5, Gh[[2, 2]]}}]

(*U系視点の表示*)
PuAb = Table[X2Up[PxAb[[ix]], Be] // N, {ix, Length[PxAb]}];
PuAe = Table[X2Up[PxAe[[ix]], Be] // N, {ix, Length[PxAe]}];
PuVe = Table[X2Up[PxVe[[ix]], Be] // N, {ix, Length[PxAe]}];
grVecVu = 
  Graphics[{Thick, RGBColor[0, 0.5, 0.1, 1], Arrow[Transpose[{PuAb, PuAb + (PuVe - PuAb)/d}]]}];
grVecAu = Graphics[{Thick, Blue, Arrow[Transpose[{PuAb, PuAb + (PuAe - PuAb)/d}]]}];
grVecBu = Graphics[{Thick, Orange, Arrow[Transpose[{PuBb, PuBe}]]}];
got = Graphics[{AbsoluteThickness[3], RGBColor[1, 0, 0, 1], dLine[X2Up, Be, opL3, 10],
    Table[Text[Style[opTu3[[ix, 1]], Large, FontFamily -> "Times"], opTu3[[ix, 2]]], {ix, Length[opTu3]}]}];
guH = grdMesh2D[X2X, Be, {{-1, 5}, {-4, 5}}, cGrnH, 1];
gx2uL = grdMesh2D[X2Up, Be, {{-1, 5}, {-4, 5}}, cBluL, 100];
Show[guH, gx2uL, got, grVecVu, grVecAu, grVecBu, PlotRange -> {{1, 3.7}, {0, 1.5}}]

(*P1の反変地縛ベクトル（青）を始点に設定し測地線上を、つまり速度ベクトルと一緒に移動させる*)
iPu = X2Up[opL3[[1]], Be] // N;
geo = True; TorB = True;

n = 5; dt = 1/n;
vxdt = (opL3[[2]] - opL3[[1]]) dt // N;
iVu = dxu[X2Upf, opL3[[1]], Be, vxdt];
pV = pVectTr2dVW[U2Xp, u, Be, iPu, {iVu, RGBColor[0, 0.5, 0.1, 1]}, intW[TorB], geo, n];

Show[guH, gx2uL, got, pV[[2]], pV[[3]], PlotRange -> {{1, 3.7}, {0, 1.5}}]

n = 100; dt = 1/n;
vxdt = (opL3[[2]] - opL3[[1]]) dt // N;
iVu = dxu[X2Upf, opL3[[1]], Be, vxdt];
pV = pVectTr2dVW[U2Xp, u, Be, iPu, {iVu, RGBColor[0, 0.5, 0.1, 1]}, intW[TorB], geo, n];
Show[guH, gx2uL, got, pV[[4]], pV[[5]], PlotRange -> {{1, 3.7}, {0, 1.5}}]

共変ベクトルの平行移動

(*P1の共変地縛ベクトル（オレンジ）を始点に設定し測地線上を、つまり速度ベクトルと一緒に移動させる*)
TorB = False;
iPu = X2Up[opL3[[1]], Be] // N;

n = 5; dt = 1/n;
vxdt = (opL3[[2]] - opL3[[1]]) dt // N;
iVu = dxu[X2Upf, opL3[[1]], Be, vxdt];
pV = pVectTr2dVW[U2Xp, u, Be, iPu, {iVu, RGBColor[0, 0.5, 0.1, 1]}, intW[TorB], geo, n];
Show[guH, gx2uL, got, pV[[2]], pV[[3]], PlotRange -> {{1, 3.7}, {0, 1.5}}]

n = 100; dt = 1/n;
vxdt = (opL3[[2]] - opL3[[1]]) dt // N;
iVu = dxu[X2Upf, opL3[[1]], Be, vxdt];
pV = pVectTr2dVW[U2Xp, u, Be, iPu, {iVu, RGBColor[0, 0.5, 0.1, 1]}, intW[TorB], geo, n];
Show[guH, gx2uL, got, pV[[4]], pV[[5]], PlotRange -> {{1, 3.7}, {0, 1.5}}]

任意の経路でのベクトルの平行移動

(*直前に実行した P1-C\[Rule]P2 直接ルートの結果を保存*)
pV4c = pV[[4]];
pV5c = pV[[5]];
Wu = d (PuBe[[1]] - PuAb[[1]]);

(*P1-A-B\[Rule]P2 遠回りルート*)
(*　A　*)
vxdt = (opL3[[3]] - opL3[[1]]) dt // N;
(*　ベクトルの表示が重ならない様に始点を少しずらす　*)
iPu = X2Up[opL3[[1]], B] + {-0.01, -0.02} // N;
iVu = dxu[X2Upf, opL3[[1]], B, vxdt];
pV = pVectTr2dVW[U2Xp, u, B, iPu, {iVu, RGBColor[0, 0.5, 0.1, 1]}, {Wu, Red, TorB}, geo, n];
pV4a = pV[[4]];
pV5a = pV[[5]];
(*最終点のベクトル位置を次の移動開始位置にセット　pV[[1,3]]=>Wu　*)
Wu = pV[[1, 3]];

(*　B　*)
vxdt = (opL3[[2]] - opL3[[3]]) dt // N;
(*　ベクトルの表示が重ならない様に始点を少しずらす　*)
iPu = X2Up[opL3[[3]], B] + {0.02, 0} // N;
iVu = dxu[X2Upf, opL3[[3]], B, vxdt];
pV = pVectTr2dVW[U2Xp, u, B,iPu, {iVu, RGBColor[0, 0.5, 0.1, 1]}, {Wu, Red, TorB}, geo, n];
pV4b = pV[[4]];
pV5b = pV[[5]];

Show[guH, gx2uL, got, got, pV4a, pV5a, pV4b, pV5b, pV4c, pV5c,
 PlotRange -> {{1, 3.7}, {0, 1.5}}

(*測地線のケース　X系で四角の二辺移動*)
geo = True;
n = 100;
dt = 1/n;
dxA = {0.5, 0.};
(*第1辺*)
vxdt = ({2, 1} - {1, 1}) dt // N;
iPu = X2Up[{1, 1}, B] // N;
iVu = dxu[X2Upf, {1, 1}, B, vxdt];
TorB = True;
Wu = dxu[X2Upf, {1, 1}, B, d {-0.5, 0.5}];
pV = pVectTr2dVW[U2Xp, u, B, 
   iPu, {iVu, RGBColor[0, 0.5, 0.1, 1]}, {Wu, Blue, TorB}, geo, n];
pV4a = pV[[4]];
pV5a = pV[[5]];
(*ベクトルの載せ替え*)
Wu = pV[[1, 3]];
vxdt = ({2, 2} - {2, 1}) dt // N;
iPu = X2Up[{2, 1} + {0.01, 0.01}, B] // N;
iVu = dxu[X2Upf, {2, 1}, B, vxdt];
pV = pVectTr2dVW[U2Xp, u, B,iPu, {iVu, RGBColor[0, 0.5, 0.1, 1]}, {Wu, Blue, TorB}, geo, n];
pV4b = pV[[4]];
pV5b = pV[[5]];

Show[guH, gx2uL, pV4a, pV5a, pV4b, pV5b,
 PlotRange -> {{1, 3.7}, {0, 1.5}}
 ]

(*第２辺*)
vxdt = ({1, 2} - {1, 1}) dt // N;
iPu = X2Up[{1, 1} + {0.01, 0.01}, B] // N;
iVu = dxu[X2Upf, {1, 1}, B, vxdt];
TorB = True;
Wu = dxu[X2Upf, {1, 1}, B, d {-0.5, 0.5}];
pV = pVectTr2dVW[U2Xp, u, B, 
   iPu, {iVu, RGBColor[0, 0.5, 0.1, 1]}, {Wu, Red, TorB}, geo, n];
pV4c = pV[[4]];
pV5c = pV[[5]];
(*ベクトルの載せ替え*)
Wu = pV[[1, 3]];
vxdt = ({2, 2} - {1, 2}) dt // N;
iPu = X2Up[{1, 2} + {0, 0.}, B] // N;
iVu = dxu[X2Upf, {1, 2}, B, vxdt];
pV = pVectTr2dVW[U2Xp, u, B,iPu, {iVu, RGBColor[0, 0.5, 0.1, 1]}, {Wu, Red, TorB}, geo, n];
pV4d = pV[[4]];
pV5d = pV[[5]];

Show[guH, gx2uL, pV4a, pV5a, pV4b, pV5b, pV4c, pV5c, pV4d, pV5d,
 PlotRange -> {{1, 3.7}, {0, 1.5}}
 ]

(*非測地線のケース　U系で四角の二辺移動*)
n = 100;
dt = 1/n;
dxA = {0.5, 0.};
(*第1辺*)
iPu = {1.5, 0.5} + {0.01, -0.01};
iVu = dt {1, 0};
TorB = True;
Wu = dxu[X2Upf, iPu, B, d {0, 0.5}];
geo = False;
pV = pVectTr2dVW[U2Xp, u, B, 
   iPu, {iVu, RGBColor[0, 0.5, 0.1, 1]}, {Wu, Blue, TorB}, geo, n];
pV4a = pV[[4]];
pV5a = pV[[5]];
(*ベクトルの載せ替え*)
Wu = pV[[1, 3]];
iPu = iPu + {1, 0} + {0.025, 0};
iVu = dt {0, 0.5};
pV = pVectTr2dVW[U2Xp, u, B,iPu, {iVu, RGBColor[0, 0.5, 0.1, 1]}, {Wu, Blue, TorB}, geo, n];
pV4b = pV[[4]];
pV5b = pV[[5]];
Show[guH, gx2uL, pV4a, pV5a, pV4b, pV5b,
 PlotRange -> {{1, 3.7}, {0, 1.5}}
 ]

(*第２辺*)
iPu = {1.5, 0.5} + {-0.01, 0};
iVu = dt {0, 0.5};
Wu = dxu[X2Upf, iPu, B, d {0, 0.5}];
pV = pVectTr2dVW[U2Xp, u, B, 
   iPu, {iVu, RGBColor[0, 0.5, 0.1, 1]}, {Wu, Red, TorB}, geo, n];
pV4c = pV[[4]];
pV5c = pV[[5]];
(*ベクトルの載せ替え*)
Wu = pV[[1, 3]];
iPu = iPu + {0, 0.5} + {0.02, 0};
iVu = dt {1, 0};
pV = pVectTr2dVW[U2Xp, u, B, 
   iPu, {iVu, RGBColor[0, 0.5, 0.1, 1]}, {Wu, Red, TorB}, geo, n];
pV4d = pV[[4]];
pV5d = pV[[5]];

Show[guH, gx2uL, pV4a, pV5a, pV4b, pV5b, pV4c, pV5c, pV4d, pV5d,
 PlotRange -> {{1, 3.7}, {0, 1.5}}
 ]

2.2 曲率

(* (+)Phi-(-)Theata Grid *)
map32[u_] := u[[{2, 1}]].{{1, 0}, {0, -1}};
X = {x1, x2, x3}; U = {u1, u2, u3};
Be = DiagonalMatrix[{6/Pi, 6/Pi}];

grRange = {{0, 4}, {-4, 0}};
scText = {{1, "/6 \[Pi]"}, {-1, "/6 \[Pi]"}};
scOfst = {{0.05, 0.1}, {-0.17, 0.0}};
axVectors = {{{0, 0}, {4, 0}}, {{0, 0}, {0, -4}}};
axTexts = {"\[Phi]", "\[Theta]"};
axOfst = {{-0.12, -0.12}, {0.15, 0.12}};
cBluH = RGBColor[0, 0.1, 0.5, 1]; cBluL = RGBColor[0, 0.1, 0.5, 0.2];
cGrnH = RGBColor[0, 0.5, 0.1, 1]; cGrnL = RGBColor[0, 0.5, 0.1, 0.3];

ax = tfCrdSys2D[U2Xsc, Be, grRange, axVectors, axTexts, scText, 
   scOfst, axOfst, cGrnH, 1];


Uc = Part[U, {2, 3}];

sz = Length[U];
QxQu = Jaco[U2Xsp[U, 0], Uc];
QxQu // MatrixForm
gxxsp = Transpose[QxQu].Diag[sz].QxQu // PowerExpand // ExpandAll // 
  Simplify
Chsp = C3s[gxxsp, Uc]
R4sp = R4[Chsp, Uc]

gSph = SphericalPlot3D[1, {theta, 0, Pi}, {phi, 0, Pi}, 
   MeshStyle -> RGBColor[0, 0.1, 0.5, 1], 
   PlotStyle -> Directive[Opacity[1], RGBColor[0, 0.1, 0.5, 0.2]], 
   Boxed -> False, Axes -> False];


pFig221[n_, m_, sf_] := Module[{},
  sn = DiagonalMatrix[{1, 1}];
  fig = 1 {{0, 0}, {Pi/12, 0}, {0, 8 Pi/12}, {Pi/12, 8 Pi/12}};
  (*fig=0.001{{0,0},{1Pi/12,0},{0,1Pi/12},{1Pi/12,1Pi/12}};*)
  du1 = fig[[3]] - fig[[1]];
  du2 = fig[[4]] - fig[[3]];
  
  pntAdr = Table[sn.fig[[ix]] + sf, {ix, Length[fig]}] // N;
  pntTxt = {"P", "a", "b", "Q"};
  g2DP00P11 = 
   Graphics[{Red, Table[{Disk[map32[pntAdr[[ix]]] // N, 0.02],
       Text[Style[pntTxt[[ix]], FontSize -> 16], 
         map32[pntAdr[[ix]]] + {0.08, 0.08}] // N   }, {ix, 
       Length[pntAdr]}   ]}];
  g3DP00P11 = 
   Graphics3D[{Red, 
     Table[{Sphere[U2Xsp[Prepend[pntAdr[[ix]], 1], Be] // N, 0.02], 
       Text[Style[pntTxt[[ix]], FontSize -> 16], 
        U2Xsp[Prepend[pntAdr[[ix]], 1] + {0.15, 0, 0}, Be] // 
         N]}, {ix, Length[pntAdr]}]}];
  RLx = Table[(U2Xsp[U, 0] /. {u1 -> 1, u2 -> Pi/2, 
       u3 -> Pi*ix/100}), {ix, 100}];
  gVR = Graphics3D[{RGBColor[0, 0.5, 0.1, 1], Thick, Dashed, Red, 
     Line[RLx]}];
  sample = Range[1, n + 2, m];
  geo = False;
  
  iPu = pntAdr[[1]];
  iVu = (pntAdr[[2]] - pntAdr[[1]])*dt;
  iWu = 0.2 {Sin[6 Pi/12], Cos[6 Pi/12]};
  (*iWu={-0.2,0}*)
  
  Du1 = {1, 0} d;
  Du2 = {0, 1} d;
  A = n {1, 0};
  Rsp = Table[
     Sum[R4sp[[i, h, j, k]] A[[h]] Du1[[j]] Du2[[k]], {h, 2}, {j, 
       2}, {k, 2}], {i, 2}] /. Thread[Uc -> iPu];
  Print["Rsp=", Rsp]
   
   
   {PuList, VuList, WuList} = 
    VectTrans[{u2, u3}, Chsp, iPu, iVu , iWu, True, geo, n];
  
  (*====1 2D ====*)
  gVu0001 = 
   Graphics[{RGBColor[0, 0.5, 0.1, 1], Thick, 
     RGBColor[0, 0.5, 0.1, 1], Dashed, Arrowheads -> Medium, 
     (*Arrow[Table[{PuList[[ix]][[2]],-PuList[[ix]][[1]]},{ix,Length[
     PuList]}]]}];*)
     Arrow[Table[map32[PuList[[ix]]], {ix, Length[PuList]}]]}];
  gWus0001 = 
   Graphics[{Blue, Thick, 
     Arrow[Table[{map32[PuList[[ix]]], 
         map32[PuList[[ix]] + WuList[[ix]]]} // N, {ix, sample}]]}];
  (*Show[ ax,g2DP00P11,gVu0001,gWus0001]*)
  (*====1 3D ====*)
  Vx = Table[
    U2Xsp[U, 0] /. u1 -> 1 /. Thread[{u2, u3} -> PuList[[ix]]], {ix, 
     Length[PuList]}];
  gVx0001 = 
   Graphics3D[{RGBColor[0, 0.1, 0.5, 1], Thick, Dashed, 
     Arrowheads -> Medium, Arrow[Vx]}];
  (*Wxe=Table[U2Xsp[U,0]/.u1\[Rule]1/.Thread[{u2,u3}\[Rule]PuList[[
  ix]]+WuList[[ix]]],{ix,n+2}];*)
  Wxe = Table[
    Vx[[ix]] + QxQu.WuList[[ix]] /. 
     Thread[U -> {1, PuList[[ix, 1]], PuList[[ix, 2]]}], {ix, n + 2}];
  Wxs = Table[{Vx[[ix]], Wxe[[ix]]}, {ix, sample}];
  gWxs0001 = Graphics3D[{Arrowheads[0.02], Blue, Thick, Arrow[Wxs]}];
  (*Show[gSph,g3DP00P11,gVx0001,gWxs0001,gVR]*)
  Pu3ds = {1, PuList[[1, 1]], PuList[[1, 2]]};
  Pu3de = {1, PuList[[Length[PuList], 1]], 
    PuList[[Length[PuList], 2]]};
  
  
  Wu = WuList[[Length[WuList]]];
  iPu = pntAdr[[2]];
  iVu = (pntAdr[[4]] - pntAdr[[2]])*dt;
  {PuList, VuList, WuList} = 
   VectTrans[{u2, u3}, Chsp, iPu, iVu , Wu, True, geo, n];
  (*====2 2D ====*)
  gVu0111 = 
   Graphics[{RGBColor[0, 0.5, 0.1, 1], Thick, 
     RGBColor[0, 0.5, 0.1, 1], Dashed, Arrowheads -> Medium, 
     Arrow[Table[map32[PuList[[ix]]], {ix, Length[PuList]}]]}];
  gWus0111 = 
   Graphics[{Blue, Thick, 
     Arrow[Table[{map32[PuList[[ix]]], 
         map32[PuList[[ix]] + WuList[[ix]]]} // N, {ix, sample}]]}];
  (*====2 3D ====*)
  Vx = Table[
    U2Xsp[U, 0] /. u1 -> 1 /. Thread[{u2, u3} -> PuList[[ix]]], {ix, 
     Length[PuList]}];
  gVx0111 = 
   Graphics3D[{RGBColor[0, 0.1, 0.5, 1], Thick, Dashed, 
     Arrowheads -> Medium, Arrow[Vx]}];
  Wxe = Table[
    Vx[[ix]] + QxQu.WuList[[ix]] /. 
     Thread[U -> {1, PuList[[ix, 1]], PuList[[ix, 2]]}], {ix, n + 2}];
  Wxs = Table[{Vx[[ix]], Wxe[[ix]]}, {ix, sample}];
  gWxs0111 = Graphics3D[{Arrowheads[0.02], Blue, Thick, Arrow[Wxs]}];
  Wua = WuList[[n + 2]] /. 
    Thread[U -> {1, PuList[[n + 2, 1]], PuList[[n + 2, 2]]}];
  Print["Wua=", Wua]
   Wxa = QxQu.WuList[[n + 2]] /. 
     Thread[U -> {1, PuList[[n + 2, 1]], PuList[[n + 2, 2]]}];
  QxQua = 
   QxQu /. Thread[U -> {1, PuList[[n + 2, 1]], PuList[[n + 2, 2]]}];
  Print["Pua=", 
    Thread[U -> {1, PuList[[n + 2, 1]], PuList[[n + 2, 2]]}]]
   
   iPu = pntAdr[[1]];
  iVu = (pntAdr[[3]] - pntAdr[[1]])*dt;
  {PuList, VuList, WuList} = 
   VectTrans[{u2, u3}, Chsp, iPu, iVu , iWu, True, geo, n];
  (*====3 2D ====*)
  gVu0010 = 
   Graphics[{RGBColor[0, 0.5, 0.1, 1], Thick, 
     RGBColor[0, 0.5, 0.1, 1], Dashed, Arrowheads -> Medium, 
     Arrow[Table[map32[PuList[[ix]]], {ix, Length[PuList]}]]}];
  gWus0010 = 
   Graphics[{Red, Thick, 
     Arrow[Table[{map32[PuList[[ix]]], 
         map32[PuList[[ix]] + WuList[[ix]]]} // N, {ix, sample}]]}];
  (*====3 3D ====*)
  Vx = Table[
    U2Xsp[U, 0] /. u1 -> 1 /. Thread[{u2, u3} -> PuList[[ix]]], {ix, 
     Length[PuList]}];
  gVx0010 = 
   Graphics3D[{RGBColor[0, 0.1, 0.5, 1], Thick, Dashed, 
     Arrowheads -> Medium, Arrow[Vx]}];
  (*Wxe=Table[U2Xsp[U,0]/.u1\[Rule]1/.Thread[{u2,u3}\[Rule]PuList[[
  ix]]+WuList[[ix]]],{ix,n+2}];*)
  Wxe = Table[
    Vx[[ix]] + QxQu.WuList[[ix]] /. 
     Thread[U -> {1, PuList[[ix, 1]], PuList[[ix, 2]]}], {ix, n + 2}];
  
  Wxs = Table[{Vx[[ix]], Wxe[[ix]]}, {ix, sample}];
  gWxs0010 = Graphics3D[{Arrowheads[0.02], Red, Thick, Arrow[Wxs]}];
  (*Show[gSph,g3DP00P11,gVx0001,gWxs0001,gVx0111,gWxs0111,gVx0010,
  gWxs0010,gVR]*)
  
  Wu = WuList[[Length[WuList]]];
  iPu = pntAdr[[3]];
  iVu = (pntAdr[[4]] - pntAdr[[3]])*dt;
  {PuList, VuList, WuList} = 
   VectTrans[{u2, u3}, Chsp, iPu, iVu , Wu, True, geo, n];
  (*====4 2D ====*)
  gVu1011 = 
   Graphics[{RGBColor[0, 0.5, 0.1, 1], Thick, 
     RGBColor[0, 0.5, 0.1, 1], Dashed, Arrowheads -> Medium, 
     (*Arrow[Table[{PuList[[ix]][[2]],-PuList[[ix]][[1]]},{ix,Length[
     PuList]}]]}];*)
     Arrow[Table[map32[PuList[[ix]]], {ix, Length[PuList]}]]}];
  gWus1011 = 
   Graphics[{Red, Thick, 
     Arrow[Table[{map32[PuList[[ix]]], 
         map32[PuList[[ix]] + WuList[[ix]]]} // N, {ix, sample}]]}];
  gru = Show[ ax, g2DP00P11, gVu0001, gWus0001, gVu0111, gWus0111, 
    gVu0010, gWus0010, gVu1011, gWus1011];
  (*====4 3D ====*)
  Vx = Table[
    U2Xsp[U, 0] /. u1 -> 1 /. Thread[{u2, u3} -> PuList[[ix]]], {ix, 
     Length[PuList]}];
  gVx1011 = 
   Graphics3D[{RGBColor[0, 0.1, 0.5, 1], Thick, Dashed, 
     Arrowheads -> Medium, Arrow[Vx]}];
  Wxe = Table[
    Vx[[ix]] + QxQu.WuList[[ix]] /. 
     Thread[U -> {1, PuList[[ix, 1]], PuList[[ix, 2]]}], {ix, n + 2}];
  Wxs = Table[{Vx[[ix]], Wxe[[ix]]}, {ix, sample}];
  gWxs1011 = Graphics3D[{Arrowheads[0.02], Red, Thick, Arrow[Wxs]}];
  Wub = WuList[[n + 2]] /. 
    Thread[U -> {1, PuList[[n + 2, 1]], PuList[[n + 2, 2]]}];
  Wxb = QxQu.WuList[[n + 2]] /. 
    Thread[U -> {1, PuList[[n + 2, 1]], PuList[[n + 2, 2]]}];
  QxQub = 
   QxQu /. Thread[U -> {1, PuList[[n + 2, 1]], PuList[[n + 2, 2]]}];
  Print["Pub=", 
   Thread[U -> {1, PuList[[n + 2, 1]], PuList[[n + 2, 2]]}]];
  
  grx = Show[gSph, g3DP00P11, gVx0001, gWxs0001, gVx0111, gWxs0111, 
    gVx0010, gWxs0010, gVx1011, gWxs1011, gVR];
  {gru, grx}]

赤道付近でのベクトルの移動経路差実験

sf = {5 Pi/12, 0};(*　移動開始点PのU座標位置 赤道近くに設定*)
n = 1000;(*　繰り返し回数 *)
dt = 1/n;
m = 500;(*　ベクトル描画サンプリング間隔 *)
res = pFig221[n, m, sf];
res[[1]]
res[[2]]

北極付近でのベクトルの移動経路差実験

sf = {0.5 Pi/12, 0};(*　移動開始点PのU座標位置 北極近くに変更*)
n = 1000; dt = 1/n;
m = 500;
res = pFig221[n, m, sf];
res[[1]]
res[[2]]

曲率テンソルRに関する実験

pFigDSvsA[U1_] := Module[{gSph, giPu, gS, gAx, gR4x},
  sz = Length[U];
  QxQu = Jaco[U2Xsp[U, 0], Uc] /. u1 -> U1;
      (*QxQu//MatrixForm;*)
  QxQu33 = Jaco[U2Xsp[U, 0], U];
  QuQx33 = Inverse[QxQu33];
  gxxsp = 
   Transpose[QxQu].Diag[sz].QxQu // PowerExpand // ExpandAll // 
    Simplify;
  Chsp = C3s[gxxsp, Uc];
  R4sp = -R4[Chsp, Uc];
  SC = 0.4;
  dd = 0.001;
  (*iPu3={u1,Pi/6,Pi/8}/.u1\[Rule]U1;*)
  iPx = U2Xsp[iPu3, 0];
  aa = {-0.5, -1};
  Au3 = {0, aa[[1]], aa[[2]]};
  dr1 = {d1, 0};
  dr2 = {0, d2};
  Dr1 = {0, d1, 0};
  Dr2 = {0, 0, d2};
  Print["QxQu33=", QxQu33 /. Thread[U -> iPu3]];
  gSph = SphericalPlot3D[U1, {u2, 0, Pi}, {u3, 0, Pi}, 
    MeshStyle -> RGBColor[0, 0.1, 0.5, 1], 
    PlotStyle -> Directive[Opacity[1], RGBColor[0, 0.1, 0.5, 0.2]], 
    Boxed -> False, Axes -> False];
  giPu = Graphics3D[{Red, Sphere[U2Xsp[iPu3, 0], 0.02]}];
  S = Norm[Cross[QxQu.dr1, QxQu.dr2]] /. d1 -> dd /. d2 -> dd /.Thread[U -> iPu3];
  Print["S=", S];
  gsc = 0.5;
  gS = Graphics3D[{Arrowheads[0.02], Black, 
        Arrow[{iPx, iPx + SC QxQu33.Dr1/dd}], Black, 
        Arrow[{iPx, iPx + SC QxQu33.Dr2/dd}],(*Green,Arrow[{iPx, iPx+
        Cross[QxQu33.Dr1/dd,QxQu33.Dr2/dd]}]*), White,
        Polygon[{iPx, iPx + SC QxQu33.Dr1/dd, 
          iPx + SC QxQu33.Dr1/dd + SC QxQu33.Dr2/dd, 
          iPx + SC QxQu33.Dr2/dd}]} /. d1 -> dd /. d2 -> dd /. 
     Thread[U -> iPu3]];
  
  Ax3 = QxQu33.Au3;
  gAx = Graphics3D[{Arrowheads[0.02], Red, Thick, 
     Arrow[{iPx, iPx + SC Ax3} /. Thread[U -> iPu3]]}];
  
  Chsp + d1 D[Chsp, u2] /. u2 -> iPu3[[2]] /. d1 -> dd // Simplify;
  Chsp /. u2 -> iPu3[[2]] + dd // Simplify;
  ChChPa = 
   Table[Sum[Chsp[[i, j, k]] dr1[[j]] aa[[k]], {j, 2}, {k, 2}], {i, 2}] /. u2 -> iPu3[[2]];
  ChChQa = aa - ChChPa - Table[Sum[Chsp[[i, j, k]] dr2[[j]] (aa - ChChPa)[[k]], {j, 2}, {k, 2}], {i, 2}]
    /. u2 -> iPu3[[2]] + d1;
  ChChPb = 
   Table[Sum[Chsp[[i, j, k]] dr2[[j]] aa[[k]], {j, 2}, {k, 2}], {i, 2}] /. u2 -> iPu3[[2]];
  ChChQb = aa - ChChPb - Table[Sum[ Chsp[[i, j, k]] dr1[[j]] (aa - ChChPb)[[k]], {j, 2}, {k, 2}], {i, 2}] 
    /. u2 -> iPu3[[2]];
  Rchch = (ChChQa - ChChQb)/S /. d1 -> dd /. d2 -> dd // Simplify;
  Print["Rchch Ra=", (ChChQa)/S /. d1 -> dd /. d2 -> dd // Simplify];
  Print["Rchch Rb=", (ChChQb)/S /. d1 -> dd /. d2 -> dd // Simplify];
  Print["Rchch=", Rchch];
  R4u = -Table[Sum[R4sp [[i, h, j, k]] (aa /. Thread[U -> iPu3])[[h]] dr1[[j]] dr2[[k]], 
    {h, 2}, {j, 2}, {k, 2}], {i, 2}]/S /. u2 -> iPu3[[2]] /. d1 -> dd /. d2 -> dd;
  Print["R4u=", R4u /. Thread[U -> iPu3]];
  R4x = QxQu33.{0, R4u[[1]], R4u[[2]]} /. Thread[U -> iPu3];
  Print["R4x=", R4x /. Thread[U -> iPu3]];
  Print["|Rx4|=", Sqrt[R4x.R4x] /. Thread[U -> iPu3]];
  Print["|Ax3|=", Sqrt[Ax3.Ax3] /. Thread[U -> iPu3]];
  Print["R4x/Ax3=", Sqrt[R4x.R4x]/Sqrt[Ax3.Ax3] /. Thread[U -> iPu3]];
  tIAI = ToString[(Sqrt[Ax3.Ax3] /. Thread[U -> iPu3]) // N];
  gR4x = Graphics3D[{Arrowheads[0.02], Red, Dashed, Thick, 
     Arrow[{iPx, iPx + SC R4x} /. Thread[U -> iPu3]], 
     Text[Style["|A|=" <> tIAI, Larger, Black, 
       FontFamily -> "Times"], (iPx + SC Ax3) /. Thread[U -> iPu3]],
     Text[Style["|Ds|=" <> ToString[Sqrt[R4x.R4x] /. Thread[U -> iPu3]], 
       Larger, Black, FontFamily -> "Times"], iPx + SC R4x],
     Text[Style["S=" <> ToString[S], Smaller, Black, FontFamily -> "Times"], 
      iPx + 0.5 SC (QxQu33.Dr1/dd + QxQu33.Dr2/dd) /. d1 -> dd /. d2 -> dd /. Thread[U -> iPu3]],
     Text[Style["|Ds|/|A|=" <> ToString[Sqrt[R4x.R4x]/Sqrt[Ax3.Ax3] /. Thread[U -> iPu3]], 
       FontSize -> 14, Red, FontFamily -> "Helvetica"], 
      iPx + 0.15 SC (QxQu33.Dr1/dd + QxQu33.Dr2/dd) /. d1 -> dd /.  d2 -> dd /. Thread[U -> iPu3]]
     }]; {gSph, giPu, gS, gAx, gR4x}
  ]

U1 = 1;
iPu3 = {u1, Pi/2, Pi/8} /. u1 -> U1;
{gSph, giPu, gS, gAx, gR4x} = pFigDSvsA[U1];
iPu3 = {u1, Pi/6, Pi/8} /. u1 -> U1;
{gSph, giPu2, gS2, gAx2, gR4x2} = pFigDSvsA[U1];

Show[gSph, giPu, gS, gAx, gR4x, giPu2, gS2, gAx2, gR4x2, 
 ViewPoint -> {Pi, Pi/2, 2}]

U1 = Sqrt[2];
iPu3 = {u1, Pi/2, Pi/8} /. u1 -> U1;
{gSph, giPu, gS, gAx, gR4x} = pFigDSvsA[U1];
iPu3 = {u1, Pi/6, Pi/8} /. u1 -> U1;
{gSph, giPu2, gS2, gAx2, gR4x2} = pFigDSvsA[U1];

Show[gSph, giPu, gS, gAx, gR4x, giPu2, gS2, gAx2, gR4x2, 
 ViewPoint -> {Pi, Pi/2, 2}]

本稿のMathematicaのコードは以下のページで使用したものです。

一般相対性理論２ [リーマン幾何学編]

Post Views: 0