diff options
| -rw-r--r-- | 03-brans-dicke-jor.mac | 267 |
1 files changed, 267 insertions, 0 deletions
diff --git a/03-brans-dicke-jor.mac b/03-brans-dicke-jor.mac new file mode 100644 index 0000000..4e3599b --- /dev/null +++ b/03-brans-dicke-jor.mac @@ -0,0 +1,267 @@ +load("setup.mac")$ +load("tensors.mac")$ +load("series.mac")$ + + +/* Equations */ + +depends(λ, [u,r,θ,ϕ]); + +/* Einstein tensor */ +for i : 1 thru 4 do for j : 1 thru 4 do ( + G_ab[i,j] : expand(ric[i,j] - ric_s*lg[i,j]/2)); + +/* Stress-energy tensor of the theory */ +for i : 1 thru 4 do for j : 1 thru 4 do ( + T_ab[i,j] : (ω/8/%pi/λ)*( + cdiff(λ,i)*cdiff(λ,j) + - lg[i,j]*(contr(ug[ρ,σ]*cdiff(λ,ρ)*cdiff(λ,σ),[ρ,σ])/2))); + +/* Einstein equations tensor */ +for i : 1 thru 4 do for j : 1 thru 4 do ( + E_ab[i,j] : expand(G_ab[i,j] - (1/λ)*( + 8*%pi*T_ab[i,j] + + cdiff(cdiff(λ, j), i) + - contr(mcs[i,j,k]*cdiff(λ,k),[k]) + - lg[i,j]*box_op(λ))))$ + +/* Einstein equations tensor with mixed indices */ +for i : 1 thru 4 do for j : 1 thru 4 do ( + uE_ab[i,j] : expand(contr(ug[i,k]*E_ab[k,j], [k])))$ + +/* solve(uE_ab[1,2]=0, diff(β,r))[1], */ + +/* Hypersurface equations */ +hs_eqs : [ +/* Determine β */ +uE_ab[1,2]=0, +/* Determine A */ +uE_ab[1,3]=0, +/* uE_ab[1,4]=0 is trivial, would determine A' otherwise */ +/* Determine B */ +uE_ab[1,1]=0 +]$ + +/* Evolution equations */ +ev_eqs : [ +/* Determine γ */ +E_ab[3,3]=0, +/* Should be fulfilled automatically (?) */ +E_ab[4,4]=0 +/* Fulfilled trivially */ +/* E_ab[3,4]=0 */ +/* E_ab[4,3]=0 */ +]$ + +/* Suplementary conditions */ +sc_eqs : [ +uE_ab[2,1]=0, +uE_ab[2,3]=0 +/* Fulfilled trivially */ +/* uE_ab[2,4]=0 */ +]$ + +hs_eqs : trigreduce(hs_eqs)$ +sc_eqs: trigreduce(sc_eqs)$ + +hs_eqs : subst([diff(λ,ϕ)=0], hs_eqs)$ +ev_eqs : subst([diff(λ,ϕ)=0], ev_eqs)$ +sc_eqs : subst([diff(λ,ϕ)=0], sc_eqs)$ + +/* Scalar wave equation */ +sw_eq : trigreduce(box_op(λ)=0)$ + + +/* Series ansatz */ + +expansion_order : 6; + +/* Sequence elements inherit dependance on u and θ */ +depends([c, H, M, N, γc, βc, Ac, Bc, λc], [u, θ])$ + +/* Main series ansatz */ +build_main_ansatz(n) := block([], [ + γ = c/r + neg_pow_series(γc, r, 2, n), + β = H + neg_pow_series(βc, r, 1, n), + A = neg_pow_series(Ac, r, 1, n), + B = neg_pow_series(Bc, r, -1, n), + λ = neg_pow_series(λc, r, 0, n)])$ + +main_ansatz : build_main_ansatz(expansion_order)$ + +funcs : [A, B, β, γ, β_expr, γ_expr, λ]; + + +/**/ + +/* Simplify ahead to speed up calculation */ +γβ_facts : [H=0, γc[2]=0]; +sim_facts : ev([ + /* From Equation 4 (1/r term) */ + diff(Φc[0],u)=0, + diff(Φc[0],u,2)=0, + /* From Equation 4 (1/r^2 term) */ + diff(Φc[0],θ)=0, + diff(Φc[0],θ, 2)=0, + diff(diff(Φc[0],θ),u)=0, + diff(diff(Φc[0],u),θ)=0], + eval,diff); + +coeffs : append(γβ_facts, sim_facts); + +/* Equation 1 : determine β */ +e_hs_eqs[1] : expand_equation(expand(exp(2*β)*λ^2*hs_eqs[1]), coeffs)$ +ce_hs_eqs[1] : coeff_equations(e_hs_eqs[1])$ +β_coeffs : block( + [vals : []], + for i : 1 thru 5 do ( + vals : append(vals, solve(subst(vals, ce_hs_eqs[1][i][2]), βc[i]))), + vals)$ + +coeffs : append(coeffs, β_coeffs)$ + +/* Equation 2 : determine A */ +e_hs_eqs[2] : expand_equation(expand(exp(2*β)*λ^2*hs_eqs[2]), coeffs)$ +/* BUG: diff(λc[0],θ)=0 does not get substituted */ +e_hs_eqs[2] : subst(coeffs, e_hs_eqs[2])$ +ce_hs_eqs[2] : coeff_equations(e_hs_eqs[2])$ + +A_coeffs : block( + [vals : []], + for i : 1 thru expansion_order do ( + if not(i=3) then ( + vals : append(vals, solve(subst(vals, ce_hs_eqs[2][i][2]), Ac[i])))), + ev(vals, eval, diff, expand))$ + +coeffs : append(coeffs, A_coeffs)$ + +/* Equation 3 : determine B */ +print_progress : true; +e_hs_eqs[3] : expand_equation(expand(exp(2*β)*λ^2*hs_eqs[3]), coeffs)$ +/* TODO: 17 31 34 slow due to squares? */ +/* BUG: diff(λc[0],θ)=0 does not get substituted; due to missing ev? */ +e_hs_eqs[3] : subst(coeffs, e_hs_eqs[3])$ +ce_hs_eqs[3] : coeff_equations(e_hs_eqs[3])$ +B_coeffs : block( + [vals : []], + for i : 1 thru expansion_order do ( + if i=2 then ( + vals : append(vals, [Bc[0]=-2*M])) + else ( + vals : append(vals, solve(subst(vals, ce_hs_eqs[3][i][2]), Bc[i-2])))), + ev(vals, eval, diff, expand))$ + +coeffs : append(coeffs, B_coeffs)$ + + +/* Evolution equations */ + +depends(D, [u, θ]); + +/* Equation 4 : determine λ */ +e_sw_eq : expand_equation(sw_eq, coeffs)$ +ce_sw_eq : coeff_equations(e_sw_eq)$ +λu_coeffs : block( + [vals : append(sim_facts, [λc[1]=D])], + for i : 3 thru 6 do ( + vals : append(vals, solve(subst(vals, ce_sw_eq[i][2]), diff(λc[i-1],u)))), + vals)$ + +coeffs : append(coeffs, λu_coeffs)$ + +depends(C, [u, θ]); +extra_coeffs : [γc[3]=C-c^3/6, Ac[3]=2*N + 3*c*diff(c,θ) + 4*c^2*cot(θ)]$ +coeffs : append(coeffs, extra_coeffs)$ + +/* Equation 5 : determine γ */ +e_ev_eqs[2] : expand_equation(expand(exp(2*β)*exp(2*γ)*λ^2*ev_eqs[2]/sin(θ)^2), coeffs)$ +/* BUG: diff(λc[0],θ)=0 does not get substituted; due to missing ev? */ +e_ev_eqs[2] : subst(coeffs, e_ev_eqs[2])$ +ce_ev_eqs[2] : coeff_equations(e_ev_eqs[2])$ + +eq_quad : trigreduce(expand(trigsimp(subst(coeffs, + 4*solve(ce_ev_eqs[2][2][2], diff(C, u))[1])))); + +/* Relation of c to c_AB */ +ce_ev_eqs[2][1]; /* TODO: D_u * c = c_u * D ??? */ + + +/* Supplementary conditions */ + +/* Supplementary Condition 1 */ +e_sc_eqs[1] : expand_equation(expand(λ^2*sc_eqs[1]), coeffs)$ +/* BUG */ +e_sc_eqs[1] : subst(coeffs, e_sc_eqs[1])$ +eq_mono : block( + [eq : coeff(e_sc_eqs[1], r, -2)=0], + trigreduce(expand(solve(eq, diff(M,u))))); + +/* Supplementary Condition 2 */ +e_sc_eqs[2] : expand_equation(sc_eqs[2], coeffs); +eq_dipo : block( + [eq : coeff(e_sc_eqs[2], r, -2)=0], + trigreduce(expand(-3*solve(eq, diff(N,u))))); + + +/* Cache */ + +path(name) := printf(false, "cache/03-brans-dicke-jor/~a.lisp", name); + +save(path("e_hs_eqs"), e_hs_eqs); +save(path("ce_hs_eqs"), ce_hs_eqs); +save(path("e_sw_eq"), e_sw_eq); +save(path("ce_sw_eq"), ce_sw_eq); +save(path("e_ev_eqs"), e_ev_eqs); +save(path("ce_ev_eqs"), ce_ev_eqs); + +save(path("beta_coeffs"), β_coeffs); +save(path("a_coeffs"), A_coeffs); +save(path("b_coeffs"), B_coeffs); +save(path("lambdau_coeffs"), λu_coeffs); + +save(path("e_sc_eqs"), e_sc_eqs); +save(path("coeffs"), coeffs) + +save(path("eq_quad"), eq_quad); +save(path("eq_mono"), eq_mono); +save(path("eq_dipo"), eq_dipo); + +load_cache() := block( + [cache : [ + "e_hs_eqs", "ce_hs_eqs", + "e_sw_eq", "ce_sw_eq", + "e_ev_eqs", "ce_ev_eqs", + "beta_coeffs", "a_coeffs", "b_coeffs", "lambdau_coeffs", + "e_sc_eqs", + "coeffs", + "eq_quad", "eq_mono", "eq_dipo"]], + for file in cache do load(path(file))); + + +/* Misc. checks */ + +p1 : (4/r + 2*diff(λ,r)/λ)*diff(β,r); +p2 : - diff(λ,r,2)/λ; +p3 : - ω*diff(λ,r)^2/λ^2; +p4 : - sum(sum(sum(sum(uh_AB[A,B]*uh_AB[C,D]*diff(h_AB[A,C],r)*diff(h_AB[B,D],r), A, 3, 4), B, 3, 4), C, 3, 4), D, 3, 4)/4; +eq1 : p1 + p2 + p3 + p4; + +expand(λ^2 * (p1 + p2 + p3 + p4)) + +expand(λ^2 * (p1)) + +expand(exp(2*β)*λ*hs_eqs[1]) + +expand(λ^2 * (p4)) + +uE_ab[1,2]=0; + +tmp : subst([diff(diff(λc[0], u), θ)=0, diff(λc[0], u, 2)=0],eq_mono) + +tmp : trigreduce(expand(trigsimp(tmp))) + +λu_coeffs[5] + +subst(diff(λu_coeffs[5], u), ev(tmp, eval)) + +/* define eth and compare with paper */ |