Refactor Brans-Dicke code

- Eliminate spurious higher orders due to previous bug in expand_term
- Use extra constant for λ_0; avoid trivial derivatives

1 files changed, 15 insertions, 57 deletions

diff --git a/03-brans-dicke-jor.mac b/03-brans-dicke-jor.mac
index 615c024..a5f0a74 100644
--- a/03-brans-dicke-jor.mac
+++ b/03-brans-dicke-jor.mac
@@ -85,7 +85,7 @@ build_main_ansatz(n) := block([], [
   β = 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)])$
+  λ = λ_0 + neg_pow_series(λc, r, 1, n)])$
 
 main_ansatz : build_main_ansatz(expansion_order)$
 
@@ -96,25 +96,15 @@ 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);
+coeffs : γβ_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 (
+  for i : 1 thru 4 do (
     vals : append(vals, solve(subst(vals, ce_hs_eqs[1][i][2]), βc[i]))),
   vals)$
 
@@ -122,8 +112,6 @@ 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(
@@ -138,13 +126,11 @@ 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 (
+  for i : 1 thru expansion_order-1 do (
     if i=2 then (
       vals : append(vals, [Bc[0]=-2*M]))
     else (
@@ -161,10 +147,11 @@ 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 : [λc[1]=D]],
+  for i : 1 thru 4 do (
+    vals : append(vals, solve(subst(vals, ce_sw_eq[i][2]), diff(λc[i+1],u)))),
   vals)$
 
 coeffs : append(coeffs, λu_coeffs)$
@@ -175,8 +162,6 @@ 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,
@@ -190,17 +175,19 @@ ce_ev_eqs[2][1]; /* TODO: D_u * c = c_u * D ??? */
 
 /* 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)))));
+  eq : solve(eq, diff(M,u))[1],
+  trigreduce(expand(eq)));
 
 /* 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)))));
+  eq : -3*solve(eq, diff(N,u))[1],
+  trigreduce(expand(eq)));
 
 
 /* Cache */
@@ -238,32 +225,3 @@ load_cache() := block(
   for file in cache do load(path(file)));
 
 load_cache();
-
-
-/* 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 */