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

リーマン幾何学で記述された、曲がった空間でのベクトルの動きや接続、曲率などの数式は難解で、その意味するところを直観的に理解することが難しいと思う。本稿に沿って、自分で処理を確認しながら実行することにより、その数式の意味を実感できるようになると思う。
サンプルプログラム実行前の準備
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』の数式検証とグラフ化に関するプログラムになります。各章節事にリンクを貼ってあるので内容を読み比べながらプログラムを実行してください。尚、内容はほぼ同じだが数式や画像の構成に多少の違いがあることをご承知ねがいたい。
2.1 共変微分
反変ベクトルの平行移動
手順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}} ]

(*第2辺*) 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}} ]

(*第2辺*) 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のコードは以下のページで使用したものです。


コメント