[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”を読み込ます。これにより以降のプログラムで使用可能となります。

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

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

prON = False;

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

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

1.2 曲がった座標系

極座標変換

手順１

X系とU系の位置ベクトル変数を定義し、ダウンロード済みの極座標変換関数 X2Up 逆変換関数 X2Upの内容を確認。

x = {x1, x2};
u = {u1, u2};
Print["x= ", x // mf, ",   u= ", u // mf]
Print["Bases \!\(\*SuperscriptBox[\(C\), \(-1\)]\)= ", U2Xp[u, Be] // mf]
?X2Up
?U2Xp

手順１の結果

手順２

ｘ系を青で描画しその上に赤で三角形を表示。

(*表示図形の設定*)
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.8, 2}},
            {"P1", {2.2, 0.7}},
            {"P2", {0.7, 3.3}}
   };
opTu3 = Table[{opT3[[ix, 1]], X2Up[opT3[[ix, 2]], B]}, {ix, 
    Length[opT3]}];
(*X系における図形の描画情報生成*)
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]}]
    }];
(*目盛り文字位置オフセット　{sc1，ｓｃ2} *)
scOfst = {-0.1 {1, 1}, -0.1 {1, 1}};

gxH = Show[grdMesh2D[X2X, Be, {4 {-1, 1}, {-1, 4}}, cBluH, 1],
   grdScText2D[X2X, Be, {4 {-1, 1}, {-1, 3}}, {{1, ""}, {1, ""}}, 
    scOfst, cBluH, n],
   grdAx2D[X2X, 
    Be, {2 {-1, 1}, 
     3 {-1, 1}}, {{{0, 0}, {1, 0}}, {{0, 0}, {0, 
       1}}}, {"\!\(\*SuperscriptBox[\(x\), \(1\)]\)", 
     "\!\(\*SuperscriptBox[\(x\), \(2\)]\)"}, {{0.2, 0.3}, {0.3, 
      0.2}}, cBluH, 1]
   ];

Show[gxH, goh]

手順２の結果

手順３

ｘ系にさらに緑のU系（極座標）を重ねて表示

gxL = Show[grdMesh2D[X2X, Be, {4 {-1, 1}, {-1, 4}}, cBluL, 1]];

(*sc1 {変換後オフセット,　変換前オフセット}*)
sc1ofst = {0.1 {1, 1}, 0 {1, 1}};
(*sc2 {変換後オフセット,　変換前オフセット}*)
sc2ofst = {0 {1, -1}, 1 {1.2, 0.2}};
scOfst = {sc1ofst, sc2ofst};

gu2xH = Show[
   grdMesh2D[U2Xp, U2Xp, {4 {0, 1}, 1 Pi {0, 1}}, cGrnH, 100],
   grdScText2D[U2Xp, Be, {{1, 4}, {1, 3}}, {{1, ""}, {1, ""}}, scOfst,
     cGrnH, 100],
   grdAx2D[U2Xp, 
    Be, {{-2, 3}, {0, 
      4}}, {{{0, 0}, {1, 0}}, {{1, 0}, {1, 
       1}}}, {"\!\(\*SuperscriptBox[\(u\), \(1\)]\)", 
     "\!\(\*SuperscriptBox[\(u\), \(2\)]\)"}, {{0.2, 0.3}, {0.3, 
      0.2}}, cGrnH, 100]];
Show[gu2xH, gxL, goh]

手順３の結果

手順４

U系（極座標）に移動し、X系の座標と図形を重ねて表示。

guH = Show[grdMesh2D[X2X, Be, {{0, 4}, {-1, 2}}, cGrnH, 1],
   grdScText2D[X2X, Be, {4 {-1, 1}, {-1, 3}}, {{1, ""}, {1, ""}}, {-0.1 {1, 1}, 0.1 {1, 1}}, cGrnH, n],
   grdAx2D[X2X, Be, {2 {-1, 1}, 3 {-1, 1}}, {{{0, 0}, {1, 0}},
　　 {{0, 0}, {0, 1}}}, {"\!\(\*SuperscriptBox[\(u\), \(1\)]\)", 
     "\!\(\*SuperscriptBox[\(u\), \(2\)]\)"}, {{0.2, 0.3}, {0.3, 0.2}}, cGrnH, 1]
   ];

(*sc1 {変換後オフセット,　変換前オフセット}*)
sc1ofst = {0 {1, 1}, 1 {-0.1, 0.8}};
(*sc2 {変換後オフセット,　変換前オフセット}*)
sc2ofst = {0 {1, -1}, 1 {1.2, 0.2}};
scOfst = {sc1ofst, sc2ofst};

gx2uL = Show[grdMesh2D[X2Up, Be, {{-1, 5}, {-4, 5}}, cBluL, 100],
   grdScText2D[X2Up, Be, {{1, 3}, {1, 3}}, {{1, ""}, {1, ""}}, scOfst, cBluH, 1000],
   grdAx2D[X2Up, Be, {{-2, 3}, {0, 4}}, {{0 {1, 1}, {1, 0}}, {{0, 0.01}, {0, 1}}},
    {"\!\(\*SuperscriptBox[\(x\), \(1\)]\)", "\!\(\*SuperscriptBox[\(x\), \(2\)]\)"}, 
    {{-0.15, 0.15}, {-0.15, -0.22}}, cBluL, 100]];
(*図形情報をX\[RightArrow]U変換しU系における図形の描画情報生成*)
got = Graphics[{AbsoluteThickness[3], RGBColor[1, 0, 0, 1],
    dLine[X2Up, B, opL3, 10],
     Table[
     Text[Style[opTu3[[ix, 1]], Large, FontFamily -> "Times"], 
      opTu3[[ix, 2]]], {ix, Length[opTu3]}]
    }];

Show[gx2uL, guH, got, PlotRange -> {{0, 4.2}, {-0.2, 1.7}}]

手順４の結果

青のX系が薄く表示されているが、図形と共にかなり歪んでいる。直交表示されたX系と見比べて位置関係を確認すると歪み方が理解できて面白い。

U系での二乗距離 s²

手順１

続いてヤコビ行列を生成、メトリックｇの生成と内容を確認

j = J[U2Xp[u, Be], u] // Simplify;
Print["J= ", j // mf];
gxx = Transpose[j].delta.j // Simplify;
Print["Metric g= \!\(\*SuperscriptBox[\(J\), \(T\)]\).\[Delta].J = ", 
Transpose[j] // mf, delta // mf, j // mf, " = ", gxx // mf]

手順１の結果

手順２

メトリックｇを使ってP1－P2の分割線分（1分割、4分割、1000分割）を引く。各線分の総計（総距離）を表示

nt = 1;
calcdS[nt_] := Module[{},
  dP = dPoint[X2Up, Be, {{2, 1}, {1, 3}}, nt][[1]];
  dMet = Table[
    gxx /. u[[1]] -> dP[[ix, 1]] /. u[[2]] -> dP[[ix, 2]], {ix, 1, 
     Length[dP] - 1}];
  dS = Table[dP[[ix + 1]] - dP[[ix]], {ix, 1, Length[dP] - 1}];
  S2 = nt*Sum[dMet[[ix]].dS[[ix]].dS[[ix]], {ix, 1, Length[dP] - 1}] //
     N;
  Graphics[{Red, Thick, Line[dP],
    Text[Style["P1", Large, FontFamily -> "Times"], 
     dP[[1]] + {0.05, -0.1}],
    Text[Style["P2", Large, FontFamily -> "Times"], 
     dP[[nt + 1]] + {0.1, 0}],
    Text[Style[
      "\!\(\*SuperscriptBox[\(s\), \(2\)]\)=" <> ToString[S2], Medium,
       Black, FontFamily -> "Times"], {3, 0.6}], 
    Text[Style["n=" <> ToString[nt], Medium, Black, 
      FontFamily -> "Times"], {2.96, 0.75}]
    }]
  ]
g1 = calcdS[1];
g2 = calcdS[4];
g3 = calcdS[1000];

Show[GraphicsGrid[{{g1, g2, g3}}], 
 Background -> RGBColor[0, 0.5, 0.1, 0.2]]

手順２の結果

1.3 ベクトルの平行移動

速度ベクトルの2次元平面移動例

(*座標描画情報を生成*)
guH = grdMesh2D[X2X, Be, {{-1, 5}, {-4, 5}}, cGrnH, 1];
gx2uL = grdMesh2D[X2Up, Be, {{-1, 5}, {-4, 5}}, cBluL, 100];

pFigPL[n_] := Module[{},
   (*ｎ分割微小時間の設定*)
   dt = 1/n;
   (*ｎ分割微小距離の計算　P1とP2の間をｎ分割*)
   vxdt = (opL3[[2]] - opL3[[1]]) dt // N;
   (*U系での初期ベクトルの始点を計算する*)
   iPu = X2Up[opL3[[1]], Be] // N;
   (*U系での初期ベクトルの終点を計算する*)
   iVu = dxu[X2Upf, opL3[[1]], Be, vxdt];
   (*pVectTr2dVWによりｎ回繰り返し平行移動を計算　ｐV[[2]]/ｐV[[2]]にｎ回分のベクトル描画情報を保存*)
   TorB = False;
   pV = pVectTr2dVW[U2Xp, u, Be, 
     iPu, {iVu, RGBColor[0, 0.5, 0.1, 1]}, {{0, 0}, Orange, TorB}, 
     True, n]
   ];
(*pFigPLで10回の平行移動を実行。　分割ベクトル描画情報 ｐV[[2]] と背景情報を描画*)
Show[gx2uL, guH, got, pFigPL[10][[2]], 
 PlotRange -> {{1, 3.7}, {0, 1.5}}]

(*pFigPLで100回の平行移動を実行。 連続ベクトル描画情報 ｐV[[4]] と背景情報を描画*)
Show[gx2uL, guH, got, pFigPL[100][[4]], PlotRange -> {{1, 3.7}, {0, 1.5}}]

速度ベクトルの球面移動例

U2Xspは基本-描画パッケージに含まれる球面変換関数を参照のこと。

(*3次元座標系の設定　U系は半径u1=２の球面{u2,u3}*)
X = {x1, x2, x3}; U = {1, u2, u3};
(*球面変換関数　U2Xsp　*)
TF = U2Xsp;

(*赤道線、P1,P2の描画情報生成*)
gRE = ParametricPlot3D[U[[1]] { Cos[p],  Sin[p], 0}, {p, 0, 2 Pi}, 
   PlotStyle -> {Red, Thick, Dashed}];
gP1P2sx = Graphics3D[{Red,
    Sphere[TF[{U[[1]], Pi/2, 0}, beta], 0.03 U[[1]]],
    Text[Style["P1", Medium, Bold], 
     TF[{U[[1]], Pi/2, 0} + {0.15 U[[1]], 0, 0}, beta]],
    Sphere[TF[{U[[1]], Pi/4, Pi/2} + {0., 0, 0}, beta], 0.03 U[[1]]],
    Text[Style["P2", Medium, Bold], 
     TF[{U[[1]], Pi/4, Pi/2} + {0.15 U[[1]], 0, 0}, beta]]
    }];
(*U系座標情報により半径:U[[1]]の球面の画像情報を生成*)
g3Do = ParametricPlot3D[
   Apply[TF, {U, 0}], {u2, 0, 1 Pi}, {u3, 0, 1 Pi}, 
   PlotStyle -> Directive[Opacity[2], RGBColor[0, 0.1, 0.5, 0.2]], 
   Boxed -> False, Axes -> False];

(*初期位置ベクトルiPu設定*)
iPu = {Pi/2, 0};
(*繰り返数ｎ　微小時間の設定　*)
n = 5; dt = 1/n;
(*初期速度ベクトルiVuの設定*)
iVu = {-Pi, Pi}/Sqrt[2] 1/2*dt;
(*pVectTr3dVで速度ベクトルの移動をシミュレーション　
ｖ：速度ベクトルの描画情報、PuListB2：速度ベクトルのU系での経路情報*)
{v, PuListB2} = Part[pVectTr3dV[U, TF, Be, iPu, iVu, n], {1, 3}];
(*全てを描画*)
Show[g3Do, v, gRE, gP1P2sx, Boxed -> False]

n = 1199; dt = 1/n;
(*初期速度ベクトルiVuの設定*)
iVu = {-Pi, Pi}/Sqrt[2] 1/2*dt;
(*pVectTr3dVで速度ベクトルの移動をシミュレーション　
ｖ：速度ベクトルの描画情報、PuListB2：速度ベクトルのU系での経路情報*)
{v, PuListB2} = Part[pVectTr3dV[U, TF, Be, iPu, iVu, n], {2, 3}];
szPuListB2 = Length[PuListB2];
(*全てを描画*)
Show[g3Do, v, gRE, gP1P2sx, Boxed -> False]

(* 球面U系　u2:Phi-u3:Theata 平面でのベクトルの移動を再生 *)

(*U系の座標情報を設定*)
grRane = {{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 = {"\[Theta]", "\[Phi]"};
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];
(*beta：U系表示グリッドのスケール補正情報*)
beta = DiagonalMatrix[{6/Pi, 6/Pi}];
(*U系の座標を生成*)
ax = tfCrdSys2D[U2Xsc, beta, grRange, axVectors, axTexts, scText, 
   scOfst, axOfst, cGrnH, 1];
gP1P2su = Graphics[{Red,
    Disk[{0, -Pi/2}, 0.04],
    Text[Style["P1", Medium, Bold], {0, -Pi/2} + {0.12, -0.12}],
    Disk[{3 Pi/6, -Pi/4}, 0.04],
    Text[Style["P2", Medium, Bold], {3 Pi/6, -Pi/4} + {0.12, -0.12}],
    Thick, Dashed, Line[{{0, -Pi/2}, {4 Pi/6, -Pi/2}}]
    }];
(*速度ベクトルの保存済み経路情報PuListB2を元に連続描画用矢印を生成*)
gV2d = Graphics[{RGBColor[0, 0.5, 0.1, 1], Thick, 
    Arrowheads -> Medium, 
    Arrow[
     Table[{PuListB2[[ix]][[2]], -PuListB2[[ix]][[1]]}, {ix, 
       szPuListB2}]]}];
Show[ ax, gP1P2su, gV2d]

(*速度ベクトルの保存済み経路情報PuListB2を元に分割描画用矢印を生成*)
nn = 200;
gV2dAr = Graphics[{RGBColor[0, 0.5, 0.1, 1], Thick, 
    Arrowheads -> Medium, 
    Arrow[Table[
      {
       {PuListB2[[ix]][[2]], -PuListB2[[ix]][[1]]},
       {PuListB2[[ix + nn]][[2]], -PuListB2[[ix + nn]][[1]]}
       },
      {ix, 1, szPuListB2 - 1, nn}]]}];
TxDs = Graphics[Text[Style["ds", Medium, Bold], {1, -0.8}]];
Show[ ax, gP1P2su, gV2dAr, TxDs]

1.4 測地線

(*経路ずれの補足線の描画情報を生成*)
pusz = szPuListB2;
bez = {{0, 0}, {0.2, 0}, {0.4, 0.3}, {0.8, 0}, {1, 0}};
dif1 = BezierFunction[bez];
PuListB2p = 
  Table[PuListB2[[ix]] + {dif1[ix/pusz][[2]], 0}, {ix, pusz}];
PuListB2m = 
  Table[PuListB2[[ix]] - {dif1[ix/pusz][[2]], 0}, {ix, pusz}];

gV2dl = Graphics[{RGBColor[0, 0.5, 0.1, 1], Thick, 
    Line[Table[{PuListB2[[ix]][[2]], -PuListB2[[ix]][[1]]}, {ix, szPuListB2}]]}];
gV2dp = Graphics[{RGBColor[0, 0.5, 0.1, 1], Thick, Dotted, 
    Line[Table[{PuListB2p[[ix]][[2]], -PuListB2p[[ix]][[1]]}, {ix, szPuListB2}]]}];
gV2dm = Graphics[{RGBColor[0, 0.5, 0.1, 1], Thick, Dotted, 
    Line[Table[{PuListB2m[[ix]][[2]], -PuListB2m[[ix]][[1]]}, {ix, szPuListB2}]]}];

gDpu = Graphics[{Thin, Arrowheads -> Small,
    Arrow[{
      {PuListB2[[Round[pusz/2]]][[2]], -PuListB2[[Round[pusz/2]]][[1]]},
      {PuListB2m[[Round[pusz/2]]][[2]], -PuListB2m[[Round[pusz/2]]][[1]]}}],
    Text[Style["\[Delta]u", Medium, Bold], {PuListB2m[[Round[pusz/2]]][[2]],
    -PuListB2m[[Round[pusz/2]]][[1]]} + {-0.03, 0.03}]
    }];

Show[ ax, gP1P2su, gV2dl, gV2dp, gV2dm, gDpu]

PuListB2pは前処理での保存済み測地線（速度ベクトルの軌跡）のリスト。PuListB2pは（＋）ベジェ変形済み、PuListB2pは　（－）ベジェ変形済みのリスト。

2次元U平面上での δu.g.δu の総計の算出と3次元X空間を描画し、その δx.δx の総計算出。

sz = szPuListB2 - 1;
Uc = Part[U, {2, 3}];
QxQu = Jaco[Apply[U2Xsp, {U, Be}], Uc];
gxx = Transpose[QxQu].Diag[3].QxQu //PowerExpand //ExpandAll //Simplify;

Print["S(theory)=  ", Pi/2 // N]

Sum[(PuListB2[[ix + 1]] - PuListB2[[ix]]).gxx.(PuListB2[[ix + 1]] - PuListB2[[ix]])
    /. Thread[Uc -> PuListB2[[ix]]] //Sqrt, {ix, sz - 1}];
Print["Sm(0\[Delta]u )=  ", %]
Sum[(PuListB2p[[ix + 1]] - PuListB2p[[ix]]).gxx.(PuListB2p[[ix + 1]] - PuListB2p[[ix]])
    /. Thread[Uc -> PuListB2p[[ix]]] //Sqrt, {ix, sz - 1}];
Print["Sm(+\[Delta]u )=  ", %]
Sum[(PuListB2m[[ix + 1]] - PuListB2m[[ix]]).gxx.(PuListB2m[[ix + 1]] - PuListB2m[[ix]])
    /. Thread[Uc -> PuListB2m[[ix]]] //Sqrt, {ix, sz - 1}];
Print["Sm(-\[Delta]u )=  ", %]

PxListB2 = Table[U2Xsp[{U[[1]], PuListB2[[ix]][[1]], PuListB2[[ix]][[2]]}, Be], {ix, sz}];
PxListB2p = Table[U2Xsp[{U[[1]], PuListB2p[[ix]][[1]], PuListB2p[[ix]][[2]]}, Be], {ix, sz}];
PxListB2m = Table[U2Xsp[{U[[1]], PuListB2m[[ix]][[1]], PuListB2m[[ix]][[2]]}, Be], {ix, sz}];

geolines = 
  Graphics3D[{Thick, RGBColor[0., 0.1, 0.5, 1], Line[PxListB2],
    Dotted, Line[PxListB2p], Line[PxListB2m]}];
Show[g3Do, geolines, gP1P2sx]

Sum[(PxListB2[[ix + 1]] - PxListB2[[ix]]).(PxListB2[[ix + 1]] - PxListB2[[ix]])//Sqrt, {ix, sz - 1}];
Print["Sx(0\[Delta]u )=  ", %]
Sum[(PxListB2p[[ix + 1]] - PxListB2p[[ix]]).(PxListB2p[[ix + 1]] - PxListB2p[[ix]])//Sqrt, {ix, sz - 1}];
Print["Sx(+\[Delta]u )=  ", %]
Sum[(PxListB2m[[ix + 1]] - PxListB2m[[ix]]).(PxListB2m[[ix + 1]] - PxListB2m[[ix]])//Sqrt, {ix, sz - 1}];
Print["Sx(-\[Delta]u )=  ", %]

2次元U平面上での δu.g.δu の総計結果

S(theory)=  1.5708
Sm(0δu )=  1.57117
Sm(+δu )=  1.58184
Sm(-\δu )=  1.58191

速度ベクトル軌跡 0δu = 1.57117が理論値 π/2=1.5708とほぼ同じことを確認。ベジェ変形+δu、ベジェ変形-δuはいずれも少し長くなっていることを確認。

ベジェ変形の軌跡の3次元X空間プロット

3次元X空間 δx.δx の総計算出結果

Sx(0\[Delta]u )=  1.571
Sx(+\[Delta]u )=  1.58166
Sx(-\[Delta]u )=  1.58178

2次元U平面上の δu.g.δu とほぼ等しい事を確認。

ベジェ変形とP1-P2距離ℓの関係をプロット

SxList = {};
Do[{
  PuListB2c = 
   Table[PuListB2[[ix]] + dn {dif1[ix/pusz][[2]], 0}, {ix, pusz}];
  Sx = Sum[(PuListB2c[[ix + 1]] - PuListB2c[[ix]]).gxx.(PuListB2c[[ix + 1]] - 
         PuListB2c[[ix]]) /. Thread[Uc -> PuListB2c[[ix]]] //Sqrt, {ix, sz - 1}];
  SxList = Append[SxList, Sx]
  }, {dn, -10, 10, 0.1}]
Plot[SxList[[Round[du*10 + 100 + 1]]], {du, -10, 10}, 
 AxesOrigin -> {0, 0}, AxesLabel -> {"\[Delta]u", "l(\[Delta]u)"}, 
 PlotRange -> {{-10, 10}, {0, 3}}]

測地線が最小の極値（dℓ/dδu=0）となっていることを確認。

トーラス面の測地線

U2Xdは基本-描画パッケージに含まれるトーラス変換関数を参照のこと。

X = {x1, x2, x3}; Ud = {2, u2, u3};
TF = U2Xd;
gRE = ParametricPlot3D[(Ud[[1]] + 1) { Cos[p],  Sin[p], 0}, {p, 0, 
    2 Pi}, PlotStyle -> {Red, Thick, Dashed}];
n = 1000; dt = 1/n;
iPu = {Pi/2, 0};
iVu = {-3 Pi/8, 48 Pi/8}*dt;
g3Dot = ParametricPlot3D[
   Apply[TF, {Ud, 0}], {u2, -Pi, 0.5 Pi}, {u3, 0, 2 Pi}, 
   PlotStyle -> RGBColor[0, 0.5, 0.5, 0.7], Boxed -> False, 
   Axes -> False];
(*g3Dot=ParametricPlot3D[Apply[TF,{Ud,0}],{u2,-Pi,0.5Pi},{u3,0,2Pi},\
PlotStyle\[Rule]RGBColor[0,0.1,0.5,0.5],Boxed\[Rule]False,Axes\[Rule]\False];
*)
Show[g3Dot, pVectTr3dV[Ud, TF, Be, iPu, iVu, n][[2]], gRE, 
 Boxed -> False]