由度规算曲率标量

上周被某人说成是“不务正业”，这次就讲讲本行吧(忽悠新手，外行可以无视，同行尽管鄙视)。

当代宇宙学的理论基石是在Robertson-Walker度规下由Einstein场方程推导出的Fridemman方程，现在的绝大部分工作都是以此为起点，而且也得到了可靠的观测支持。从给定度规得出场方程的具体形式是广义相对论的基础内容，但是由于方法繁复，教科书中都不会给出具体的计算过程，而结果又不是一望便知的，学到这里谁都少不了课下的一番推导验证。我当年偷懒跳过，现在却发现自己怎么都算不对了……

怅惘之际在论文库里发现2000年《上海天文台年刊》第21期中有一篇《利用Mathematica软件表示真空Einstein场方程》的文章，但又不想为此学门新语言，便用Matlab仿写了一个，代码如下：

syms t k r theta phi;
% 声明符号变量
R=sym('R(t)');
% 声明尺度因子R为时间t的函数
g=[-1,R^2/(1-k*r^2),r^2*R^2,r^2*R^2*sin(theta)^2];
% 输入R-W度规对角元
v=[t,r,theta,phi]; 
% 输入自变量矩阵
g_mn=diag(g,0);  
%将对角元转换成4X4 矩阵
Pg=cat(3,diff(g_mn,v(1)),diff(g_mn,v(2)),diff(g_mn,v(3)),diff(g_mn,v(4)));  
% 微分度规
for i=1:4
 gama(:,:,i)=(1/g(i))/2*(squeeze(Pg(i,:,:))+squeeze(Pg(i,:,:)).'-Pg(:,:,i));  
%求Christoff联络
%虽然Pg是4X4X4的三维方阵，但Matlab中Pg(1,:,:)与Pg(:,:,1)并不相同，因此要用squeeze降维
end
G1=0;G2=0;G3=0;G4=0;  
%下面计算Ricci张量
for i=1:4
 G1=G1+jacobian(gama(i,:,i),v);
 G2=G2+diff(gama(:,:,i),v(i));
  for j=1:4
    G3=G3+gama(:,j,i)*gama(:,i,j).';
    G4=G4+gama(:,:,i)*gama(i,j,j);
  end
end
R_mn=G1-G2+G3-G4; 
[R,combine]=simple(sum(diag(R_mn)./diag(g_mn)));
%得出曲率标量R

syms t k r theta phi;

% 声明符号变量

R=sym('R(t)');

% 声明尺度因子R为时间t的函数

g=[-1,R^2/(1-k*r^2),r^2*R^2,r^2*R^2*sin(theta)^2];

% 输入R-W度规对角元

v=[t,r,theta,phi];

% 输入自变量矩阵

g_mn=diag(g,0);

%将对角元转换成4X4 矩阵

Pg=cat(3,diff(g_mn,v(1)),diff(g_mn,v(2)),diff(g_mn,v(3)),diff(g_mn,v(4)));

% 微分度规

for i=1:4

gama(:,:,i)=(1/g(i))/2*(squeeze(Pg(i,:,:))+squeeze(Pg(i,:,:)).'-Pg(:,:,i));

%求Christoff联络

%虽然Pg是4X4X4的三维方阵，但Matlab中Pg(1,:,:)与Pg(:,:,1)并不相同，因此要用squeeze降维

end

G1=0;G2=0;G3=0;G4=0;

%下面计算Ricci张量

for i=1:4

G1=G1+jacobian(gama(i,:,i),v);

G2=G2+diff(gama(:,:,i),v(i));

for j=1:4

G3=G3+gama(:,j,i)*gama(:,i,j).';

G4=G4+gama(:,:,i)*gama(i,j,j);

end

R_mn=G1-G2+G3-G4;

[R,combine]=simple(sum(diag(R_mn)./diag(g_mn)));

%得出曲率标量R

虽然编写的时候问题多多，但总算是得出正确结果了，看来Matlab的符号计算功能也不差么，为什么就没有人用呢？兴奋之余又试了一下传说中可以改善显示效果的pretty函数，看到的便是左边的结果……

我半分钟之后才反应过来这居然就是3R”/R（显然是BBS泡的太少，对ASCII码不够敏感），那一刻我终于明白为什么没人用Matlab做符号计算了……

当下决定学Mathematica，周建峰等人的代码也贴上来备查吧，shift+回车直接执行。

附：

v = {t, r, e, phi};
L1 = {-1, R[t]^2/(1 - k r^2), (r^2)(R[t]^2), (r^2)(R[t]^2)(Sin[e]^2)};
g = DiagonalMatrix[L1];
Result ={};
invg = Inverse [g];
dg1 = Outer [D,g,v];
dg2 = Transpose [dg1,{1,3,2}];
dg3 = Transpose [dg1,{2,3,1}];
G=(1/2)invg.(dg1+dg2-dg3);
dG = Outer[D,G,v];
Ruv1 = Table[Sum[dG[[n,i,n,j]],{n,4}],{i,4},{j,4}];
Ruv2 = Table[Sum[dG[[n,i,j,n]],{n,4}],{i,4},{j,4}];
Ruv3 = Table[Sum[G[[n,i,j]]G[[m,n,m]],{n,4},{m,4}],{i,4},{j,4}];
Ruv4 = Table[Sum[G[[n,i,m]]G[[m,j,n]],{n,4},{m,4}],{i,4},{j,4}];
Ruv= Ruv1-Ruv2-Ruv3+Ruv4;
Result = Append[Result,Simplify[Ruv]];
Rrr=Sum[invg[[i,i]]Ruv[[i,i]],{i,4}];
Simplify[Rrr]

v = {t, r, e, phi};

L1 = {-1, R[t]^2/(1 - k r^2), (r^2)(R[t]^2), (r^2)(R[t]^2)(Sin[e]^2)};

g = DiagonalMatrix[L1];

Result ={};

invg = Inverse [g];

dg1 = Outer [D,g,v];

dg2 = Transpose [dg1,{1,3,2}];

dg3 = Transpose [dg1,{2,3,1}];

G=(1/2)invg.(dg1+dg2-dg3);

dG = Outer[D,G,v];

Ruv1 = Table[Sum[dG[[n,i,n,j]],{n,4}],{i,4},{j,4}];

Ruv2 = Table[Sum[dG[[n,i,j,n]],{n,4}],{i,4},{j,4}];

Ruv3 = Table[Sum[G[[n,i,j]]G[[m,n,m]],{n,4},{m,4}],{i,4},{j,4}];

Ruv4 = Table[Sum[G[[n,i,m]]G[[m,j,n]],{n,4},{m,4}],{i,4},{j,4}];

Ruv= Ruv1-Ruv2-Ruv3+Ruv4;

Result = Append[Result,Simplify[Ruv]];

Rrr=Sum[invg[[i,i]]Ruv[[i,i]],{i,4}];

Simplify[Rrr]

订阅评论

4 评论

最旧