plotdf.lisp: change the name of one of the functions SUBXY which were duplicated.
[maxima.git] / share / tensor / einhil.dem
blobe6b034c50e5def6722dbc63a74c070b8ddd4fb1e
1 /* Copyright (C) 2008 Viktor T. Toth <http://www.vttoth.com/>
2  *
3  * This program is free software; you can redistribute it and/or
4  * modify it under the terms of the GNU General Public License as
5  * published by the Free Software Foundation; either version 2 of
6  * the License, or (at your option) any later version.
7  *
8  * This program is distributed in the hope that it will be
9  * useful, but WITHOUT ANY WARRANTY; without even the implied
10  * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11  * PURPOSE.  See the GNU General Public License for more details.
12  *
13  * Deriving the Einstein field equations from the Einstein-Hilbert action
14  *
15  */
16 ("Deriving the Einstein field equations in FLRW cosmology" )$
17 if get('ctensor,'version)=false then load(ctensor);
18 if get('itensor,'version)=false then load(itensor);
19 ("The first step is to construct a symmetrized Riemann tensor.")$
20 ("For this, we employ an auxiliary symmetrized metric.")$
21 remsym(g,2,0);
22 remsym(g,0,2);
23 remsym(gg,2,0);
24 remsym(gg,0,2);
25 remcomps(gg);
26 imetric(gg);
27 icurvature([a,b,c],[e])*gg([d,e],[])$
28 contract(rename(expand(%)))$
29 %,ichr2$
30 contract(rename(expand(%)))$
31 canform(%)$
32 contract(rename(expand(%)))$
33 ("We now reexpress gg by symmetrizing g using kdels and simplify:")$
34 components(gg([a,b],[]),kdels([a,b],[u,v])*g([u,v],[])/2);
35 components(gg([],[a,b]),kdels([u,v],[a,b])*g([],[u,v])/2);
36 %th(5),gg$
37 ("Some of the following simplifications may take some time...")$
38 contract(rename(expand(%th(2))))$
39 contract(canform(%))$
40 ("Now we can switch to the real metric:")$
41 imetric(g);
42 contract(rename(expand(%th(3))))$
43 ("At last, we got the covariant Riemann tensor.")$
44 remcomps(R);
45 components(R([a,b,c,d],[]),%th(3));
46 ("What we really need, though, is the curvature scalar:")$
47 g([],[a,b])*R([a,b,c,d])*g([],[c,d])$
48 contract(rename(canform(%)))$
49 contract(rename(canform(%)))$
50 components(R([],[]),%);
51 ("Before going further, we establish the symmetry properties of g.")$
52 decsym(g,2,0,[sym(all)],[]);
53 decsym(g,0,2,[],[sym(all)]);
54 ("Now we can construct the Einstein-Hilbert action.")$
55 ishow(1/(16*%pi*G)*((2*L+'R([],[])))*sqrt(-determinant(g)))$
56 L0:%,R$
57 ("We construct and simplify the 2nd order Euler-Lagrange equation:")$
58 canform(contract(canform(rename(contract(expand(diff(L0,g([],[m,n]))-idiff(diff(L0,g([],[m,n],k)),k)+idiff(rename(idiff(contract(diff(L0,g([],[m,n],k,l))),k),1000),l)))))))$
59 ishow(e([m,n],[])=canform(%*16*%pi/sqrt(-determinant(g))))$
60 ("We build a ctensor program to calculate tensor components:")$
61 EQ:ic_convert(%th(2))$
62 ("Now we set up the FLRW metric of cosmology:")$
63 ct_coords:[t,r,u,v];
64 lg:ident(4);
65 lg[2,2]:-a^2/(1-k*r^2);
66 lg[3,3]:-a^2*r^2;
67 lg[4,4]:-a^2*r^2*sin(u)^2;
68 dependencies(a(t));
69 cmetric();
70 derivabbrev:true;
71 christof(false);
72 ("We can at last evaluate the Euler-Lagrange equation...")$
73 e:zeromatrix(4,4);
74 ("This is definitely going to take a little time...")$
75 ev(EQ);
76 ("...and obtain the Einstein tensor for cosmology.")$
77 expand(radcan(ug.e));
78 ("What about the spherically symmetric case?")$
79 lg:ident(4);
80 lg[1,1]:B;
81 lg[2,2]:-A;
82 lg[3,3]:-r^2;
83 lg[4,4]:-r^2*sin(u)^2;
84 kill(dependencies);
85 dependencies(A(r),B(r));
86 cmetric();
87 christof(false);
88 e:zeromatrix(4,4);
89 ("We re-evaluate the Euler-Lagrange equation. Be patient...")$
90 ev(EQ);
91 E:expand(radcan(ug.e));
92 ("In a vacuum solution, E should be zero. Let's solve for it!")$
93 exp:findde(E,2);
94 solve(ode2(exp[1],A,r),A);
95 %,%c=-2*M;
96 a:%[1],%c=-2*M;
97 ode2(ev(exp[2],a),B,r);
98 b:ev(%,%c=rhs(solve(rhs(%)*rhs(a)=1,%c)[1]));
99 ("The solution must be consistent with the third equation")$
100 factor(ev(ev(exp[3],a,b),diff));
101 ("Finally, we obtain the Schwarzschild-de Sitter metric")$
102 expand(ev(lg,a,b));
103 /* End of demo -- comment line needed by MAXIMA to resume demo menu */