solve: do not call MEVAL.
[maxima.git] / share / tensor / bradic.dem
blobae5e8cbddc2785e84be38d65e255531669b9f41e
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 field equations from the action for Brans-Dicke cosmology
14  *
15  */
16 ("Deriving the field equations in FLRW cosmology for the Brans-Dicke action" )$
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 Brans-Dicke action.")$
55 L1:1/(16*%pi)*(('R([],[])))*sqrt(-determinant(g))$
56 L2:1/(16*%pi)*(w*gg([],[a,b])*f([],[],a)*f([],[],b)/f([],[]))*sqrt(-determinant(g))$
57 L0:radcan(-f([],[])*L1+L2)$
58 L0:contract(ev(L0,R))$
59 ("We construct and simplify the 2nd order Euler-Lagrange equation:")$
60 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)))))))$
61 ishow(e([m,n],[])=canform(%*16*%pi/sqrt(-determinant(g))))$
62 ("We build a ctensor program to calculate tensor components:")$
63 EQ:ic_convert(%th(2))$
64 ("Now we set up the FLRW metric of cosmology:")$
65 ct_coords:[t,r,u,v];
66 lg:ident(4);
67 lg[2,2]:-a^2/(1-k*r^2);
68 lg[3,3]:-a^2*r^2;
69 lg[4,4]:-a^2*r^2*sin(u)^2;
70 ("Let's not forget that the scalar field is also a function of time.")$
71 dependencies(a(t),f(t));
72 cmetric();
73 derivabbrev:true;
74 christof(false);
75 ("We can at last evaluate the Euler-Lagrange equation...")$
76 e:zeromatrix(4,4);
77 ("This is definitely going to take a little time...")$
78 ev(EQ);
79 ug.e$
80 subst(f,f([],[]),%)$
81 expand(radcan(%/f));
82 ("As a last step, let us calculate the Hamiltonian associated with f.")$
83 ic_convert(L([],[])=L0)$
84 ev(%)$
85 L:factor(subst(f,f([],[]),L));
86 P:subst(diff(f,t),ft,diff(subst(ft,diff(f,t),L),ft));
87 H:ratsimp(P*diff(f,t)-L);
88 /* End of demo -- comment line needed by MAXIMA to resume demo menu */