日志
计算黑洞度规的代码
热度 1 |||
1915年，爱因斯坦经过10年的艰苦努力，终于推导出了广义相对论的方程，又称爱因斯坦方程。不过，这个（或者这组）方程相当繁复，爱因斯坦只就一些情况作了近似解。这时正是一战鏖战激烈的时候，在德军的俄国前线，有一名炮兵中尉读到了爱因斯坦的最新论文，在炮火纷飞之中，他竟然得出了爱因斯坦方程的一个精确解。这名德国军官名叫 Karl Schwarzschild，其姓 Schwarzschild 是黑色符号的意思。
读者莫要以为德军素质奇高，其实这位炮兵原来是一名物理教授。一战那个年代，高级知识分子是没有特殊待遇的，很多杰出诗人、科学家都应征入伍，有的战死沙场，有的是战场还没看到就染病身亡。这位 Schwarzschild 第二年也病死了。
爱因斯坦方程左边是与空间时间弯曲相关的一个量，右边是空动量能量。对于真空来说，动量能量为零，所以爱因斯坦方程左边也为零，但这并不意味着空时就没有弯曲。例如，太阳之外假如不考虑其他物质，就是真空，但这些地方也是有引力的。Schwarzschild 的解就是一个真空爱因斯坦方程的解。他这个推导费了相当的劳动。
在《广义相对论计算代码》一文中，我写了几段代码从度规计算黎曼张量等。其中一个叫做 Ricci 张量。对于真空来说，这个RICCI张量为零。上次我代入 Schwarzschild 的解验算，发现 Ricci 张量确实为零。今天我想，能否用计算机解出 Schwarzschild 的结果呢？
因为球对称，我们度规是：sm = DiagonalMatrix[{ A[r], B[r], r^2, r^2 Sin[\[Theta]]^2}]
其中A、B为 r 的未知函数。
方程是 Ricci 张量为零。这是一个微分方程组
ricci = RicciT[{sm, xx}]
DSolve [ {ricci[[1, 1]] == 0, ricci[[2, 2]] == 0}, {A, B}, {r}]
结果得到：
不知怎样强迫 Mathematica 简化上面的指数与对数，我只好手动了。
A(r) = c2 - c1/r
B(r) = r c3/ (c1-r c2)
r 为无穷大时，A = -1, B =1 ，因此 c2 =1, c3 = 1
所以，
A (r) = -1 - c1/r
B (r) = 1 / ( c1/r -1）
c1 进一步由牛顿万有引力近似确定。这正是 Schwarzschild 的解。
当年学相对论有 mathematica 就好了。
/////全部代码如下
ChristoffelSym[z_] := Module[{g, x, d, ginv, res}, {g, x} = z;
d = Length[x]; ginv = Simplify[Inverse[g]];
res = Table[(1/2)*
Sum[ginv[[a,
b]]*(D[g[[b, u]], x[[v]]] + D[g[[b, v]], x[[u]]] -
D[g[[u, v]], x[[b]]]), {b, 1, d}], {a, 1, d}, {u, 1, d}, {v,
1, d}];
FullSimplify[res]]
RiemannCurvature[z_] := Module[{g, x, d, C, res}, {g, x} = z;
d = Length[x]; C = ChristoffelSym[z];
res = Table[
D[C[[a, b, v]], x[[u]]] - D[C[[a, b, u]], x[[v]]] +
Sum[C[[a, s, u]]*C[[s, b, v]], {s, 1, d}] -
Sum[C[[a, s, v]]*C[[s, b, u]], {s, 1, d}], {a, 1, d}, {b, 1,
d}, {u, 1, d}, {v, 1, d}];
Simplify[res]]
ContractMi[R_] := Module[{d, res}, d = Dimensions[R, 1][[1]];
res = Table[Sum[R[[u, a, u, b]], {u, 1, d}], {a, 1, d}, {b, 1, d}];
Simplify[res]]
KretschmannScalar[z_] := Module[{go, x, n, R, gi, res}, {go, x} = z;
R = RiemannCurvature[z];
n = Length[x];
gi = Inverse[go];
res = Sum[
R[[a, b, c, d]]*R[[e, f, g, h]]*go[[e, a]]*gi[[f, b]]*gi[[g, c]]*
gi[[h, d]], {e, 1, n}, {f, 1, n}, {g, 1, n}, {h, 1, n}, {a, 1,
n}, {b, 1, n}, {c, 1, n}, {d, 1, n}];
Simplify[res]]
RicciT[z_] := ContractMi[RiemannCurvature[z]]
RicciS2[g_, rt_] := Module[{d, ginv, res}, d = Dimensions[g, 1][[1]];
ginv = Inverse[g];
res = Sum[ginv[[u, v]]*rt[[u, v]], {u, 1, d}, {v, 1, d}];
Simplify[res]]
RicciS[z_] := Module[{x, g, d, rt, ginv, res}, {g, x} = z;
d = Length[x]; rt = RicciT[z];
ginv = Inverse[g];
res = Sum[ginv[[u, v]]*rt[[u, v]], {u, 1, d}, {v, 1, d}];
FullSimplify[res]]
EisteinTensor[z_] := Module[{x, g, d, ginv, rt, rs, res}, {g, x} = z;
d = Length[x]; ginv = Inverse[g];
rt = RicciT[z];
rs = RicciS2[g, rt];
res = rt - (1/2) g*rs;
Simplify[res]]
sm = DiagonalMatrix[{ A[r], B[r], r^2, r^2 Sin[\[Theta]]^2}]
xx = {t, r, \[Theta], \[Phi]}
ricci = RicciT[{sm, xx}]
DSolve [ {ricci[[1, 1]] == 0, ricci[[2, 2]] == 0}, {A, B}, {r}]
{A[r], B[r] } /. Out[12][[1]] // Simplify
