yuedongxiao的个人博客分享 http://blog.sciencenet.cn/u/yuedongxiao

博文

程序解爱因斯坦方程

已有 5459 次阅读 2016-4-28 15:17 |个人分类:物理|系统分类:科普集锦

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}]
结果得到:
sc.png
不知怎样强迫 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}]




https://blog.sciencenet.cn/blog-684007-973382.html

上一篇:GPS与相对论之迷途知返
下一篇:经典哲学都是兜圈子
收藏 IP: 24.7.123.*| 热度|

7 刘全慧 张云 赵建民 李颖业 张江敏 刘然 lrx

该博文允许注册用户评论 请点击登录 评论 (4 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-4-27 23:35

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部