diff options
| -rw-r--r-- | 01-bondi.mac | 165 | ||||
| -rw-r--r-- | series.mac | 296 | ||||
| -rw-r--r-- | setup.mac | 105 | ||||
| -rw-r--r-- | tensors.mac | 21 |
4 files changed, 587 insertions, 0 deletions
diff --git a/01-bondi.mac b/01-bondi.mac new file mode 100644 index 0000000..dda6bf8 --- /dev/null +++ b/01-bondi.mac @@ -0,0 +1,165 @@ +load("setup.mac")$ +load("tensors.mac")$ +load("series.mac")$ + + +/* Equations */ + +/* Stress-energy tensor of the theory */ +T_ab : zeromatrix(4,4)$ + +/* Einstein equations tensor */ +for i : 1 thru 4 do for j : 1 thru 4 do ( + E_ab[i,j] : expand(ric[i,j] - ric_s*lg[i,j]/2 - 8*%pi*T_ab[i,j]))$ + +/* 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); + + +/* Series ansatz */ + +expansion_order : 6; + +/* Sequence elements inherit dependance on u and θ */ +depends([c, H, M, N, γc, βc, Ac, Bc], [u, θ])$ + +/* Simplify ahead to speed up calculation */ +/* γβ_facts : [H=0, γc[2]=0]$ */ + +/* 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)])$ + +main_ansatz : build_main_ansatz(expansion_order)$ + +funcs : [A, B, β, γ, β_expr, γ_expr]; + + +/**/ + +/* Simplify ahead to speed up calculation */ +γβ_facts : [H=0, γc[2]=0]; + +coeffs : γβ_facts; + +/* Equation 1 : determine β */ +e_hs_eqs[1] : expand_equation(solve(hs_eqs[1], diff(β, r))[1], coeffs); +ce_hs_eqs[1] : coeff_equations(e_hs_eqs[1]); +β_coeffs : block([len, lst], + len : length(ce_hs_eqs[1]), + lst : makelist(solve(ce_hs_eqs[1][i][2], βc[i]), i, 1, len), + apply(append, lst)); + +coeffs : append(coeffs, β_coeffs); + +/* Equation 2 : determine A */ +e_hs_eqs[2] : expand_equation(expand(exp(2*β)*hs_eqs[2]), coeffs); +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 */ +e_hs_eqs[3] : expand_equation(hs_eqs[3], coeffs); +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)$ + +/* Equation 4 : determine γ (evolution equation) */ + +/* +We wish M,N,C to be related to monopole, dipole and quadrupole moments +respectively, also when converted from the static case. Thus, Ac[3] as to be of +following form. +*/ +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); + +e_ev_eqs[2] : expand_equation(ev_eqs[2], coeffs); +ce_ev_eqs[2] : coeff_equations(e_ev_eqs[2]); + +eq_quad : expand(4*solve(ce_ev_eqs[2][1][2], diff(C,u))[1]); + +/* Convoluted way to show that the results for C_u are consistent */ +e_ev_eqs[1] : expand_equation(ev_eqs[1], coeffs); +ce_ev_eqs[1] : coeff_equations(e_ev_eqs[1]); + +eq_quad_alt : trigreduce(expand(trigsimp(4*solve(ce_ev_eqs[2][2][2], diff(C,u))[1]))); +gammac4 : trigreduce(expand(trigsimp(trigreduce(expand(solve(ce_ev_eqs[1][2][2], diff(γc[4],u))))[1]))); +eq_quad_alt_res : expand(subst(gammac4, eq_quad_alt)) + + +/* Supplementary conditions */ + +/* Supplementary Condition 1 */ +e_sc_eqs[1] : expand_equation(sc_eqs[1], coeffs); +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))))); + + + +/* BUG: should not generate terms with r^-3 and so on ... */ +block([expansion_order : 2], + expand_equation(sc_eqs[1], coeffs)) + +block([expansion_order : 2], + expand_equation(sc_eqs[2], coeffs)) diff --git a/series.mac b/series.mac new file mode 100644 index 0000000..5617a1e --- /dev/null +++ b/series.mac @@ -0,0 +1,296 @@ +/* Term generation */ +exp_taylor(expr, n) := expand(subst([ξ=expr], taylor(exp(ξ), ξ, 0, n))); +neg_pow_series(c, r, a, b) := sum(c[i] / r^i, i, a, b); + + +/* Term filtering */ + +filter_order(expr, x, a, b) := block([acc, len], + acc : 0, + len : if atom(expr) or not(op(expr) = "+") then 1 else length(expr), + + for i thru len do block([deg, p], + p : if len=1 then expr else part(expr, i), + deg : hipow(partition([p], x)[2], x), + if (deg >= a and deg <= b) then acc : acc + p), + + acc); + +/* Speed up by expanding and mapping only once over the list (check operator is +) */ +keep_upto_order(expr, var, start, end) := block([parts], + expr : expand(expr), + parts : makelist(coeff(expr, var, i) * r^i, i, start, end), + apply("+", parts)); + + +/* Term generation with filtering */ + +/* +Calculate the exponential series expansion evaluated at `expr' up to `n'th order +in `x' +*/ +exp_series(expr, x, n, a, b) := block([acc, xn, len], + expr : expand(expr), + xn : 1, + acc : 1, + len : if atom(expr) or not(is(op(expr)="+")) then 1 else length(expr), + /* Assume a, b < 0 for now */ + for i thru n do block([subexpr], + /* Discard terms that would end up creating terms outside range */ + /* + subexpr : 0, + for j thru len do block([term, deg], + term : if len=1 then expr else part(expr, j), + deg : hipow(term, x), + if not((deg + lopow(xn, x)) > b or (deg + hipow(xn, x)) < a) + then + subexpr : subexpr + term), + */ + /**/ + xn : filter_order(ratexpand(xn * expr / i), x, a, b), + acc : acc + xn), + acc); + +horner_exp_series(expr, x, n, a, b) := block([acc], + /* expand(1 + x*(1 + x*(1 + x*(1+x/4)/3)/2)) */ + expr : expand(expr), + acc : 1, + for i thru n do acc : keep_upto_order(1 + expr*acc/(n-i+1), x, a, b), + acc); + + +/* Multiply series expansions upto nth order */ + +matchdeclare(m_freeof_r, lambda([e], freeof(r, e)))$ +matchdeclare(m_small, lambda([e], e < -expansion_order))$ +defrule(r_simp_small, m_freeof_r*r^m_small, 0)$ + +__simple_filter(series, max_order) := block([ret, tmp], + tmp : expansion_order, + expansion_order : max_order, + ret : apply1(series, r_simp_small), + expansion_order : tmp, + ret); + +__table_size(factors) := block( + [nfactors, ranges, t_min, t_max], + nfactors : length(factors), + ranges : make_array(any, nfactors, 2), + /* Fill table with maximal and minimal powers */ + for i : 0 thru nfactors-1 do block( + ranges[i, 0] : lopow(factors[i+1], r), + ranges[i, 1] : hipow(factors[i+1], r)), + /* Find greatest and smallest r powers */ + t_min : ranges[0, 0], + t_max : ranges[0, 1], + for i : 0 thru nfactors-1 do ( + t_min : min(t_min, ranges[i, 0]), + t_max : max(t_max, ranges[i, 1])), + [t_min, t_max]); + +__fill_term_table(table, factors, t_min, t_max, len) := block( + [nfactors], + nfactors : length(factors), + /* Split factors into term table */ + for i : 0 thru nfactors-1 do ( + /* Initialize table row with 0s */ + for j : 0 thru len-1 do table[i, j] : 0, + /**/ + if is(op(factors[i+1])="+") then ( + for term in args(factors[i+1]) do block( + [j, n], + n : lopow(term, r), + j : n - t_min, + table[i, j] : table[i, j] + coeff(term, r, n))) + else block( + [j, n, term], + term : factors[i+1], + n : lopow(term, r), + j : n - t_min, + table[i, j] : coeff(term, r, n)))); + +__mul_series_level(table, max_order, i, stack, nfactors, len, t_min) := block( + [], + if is(i = nfactors) then block( + [order, acc], + order : 0, + for j : 0 thru nfactors-1 do order : order + stack[j] + t_min, + if order >= -max_order then ( + acc : 1, + for j : 0 thru nfactors-1 do acc : acc * table[j, stack[j]], + acc * r^order) + else 0) + else block([acc], + /* TODO: omit external 0s */ + acc : 0, + for j : 0 thru len-1 do ( + stack[i] : j, + acc : acc + __mul_series_level( + table, max_order, i+1, stack, nfactors, len, t_min)), + acc)); + +/* +Multiply list of factors, each expressed as terms of series, only upto max_order +*/ +mul_series(factors, max_order) := block( + [term_table, len, t_min, t_max, stack, nfactors], + nfactors : length(factors), + + /* If there is nothing to multiple, still filter higher terms */ + if is(nfactors=1) then ( + return(__simple_filter(factors[1], max_order))), + + /* Allocate stack and table for ranges */ + stack : make_array(any, nfactors), + + /* Compute table size */ + block([pair], + pair : __table_size(factors), + t_min : pair[1], + t_max : pair[2]), + len : abs(t_max - t_min) + 1, + + /* Allocate table for terms */ + term_table : make_array(any, nfactors, len), + __fill_term_table(term_table, factors, t_min, t_max, len), + + __mul_series_level(term_table, max_order, 0, stack, nfactors, len, t_min))$ + + +/**/ + +β_expr_val(x, n) := block( + if x < 0 then ( + exp(x*H)*horner_exp_series(abs(x)*(H-rhs(main_ansatz[2])), r, n, -n, 2)) + else ( + exp(x*H)*horner_exp_series(x*(rhs(main_ansatz[2])-H), r, n, -n, 2))); + +γ_expr_val(x, n) := horner_exp_series(x*rhs(main_ansatz[1]), r, n, -n, 2); + +matchdeclare(m_any1, lambda ([e], freeof(β, γ, e))); +matchdeclare(m_any2, lambda ([e], freeof(β, γ, e))); +matchdeclare(m_any3, lambda ([e], freeof(β, γ, e))); +matchdeclare(m_beta, lambda ([e], not freeof(β, e) and freeof(γ, e))); +matchdeclare(m_gamma, lambda ([e], not freeof(γ, e) and freeof(β, e))); + +/* Replace β and γ expressions into a symbolic representation */ +defrule(r_wrap_exp_β, exp(m_any1*m_beta), β_expr(m_any1)); +defrule(r_wrap_exp_γ, exp(m_any1*m_gamma), γ_expr(m_any1)); +defrule(r_wrap_exp_βγ, exp(m_any1*m_beta + m_any2*m_gamma), + β_expr(m_any1)*γ_expr(m_any2)); +defrule(r_wrap_exp_βγ2, exp(m_any1*m_beta + m_any2*m_gamma + m_any3), + exp(m_any3)*β_expr(m_any1)*γ_expr(m_any2)); + +/* Apply wrap rules to expression */ +wrap_exp(expr) := ( + apply1(expr, r_wrap_exp_β, r_wrap_exp_γ, r_wrap_exp_βγ, r_wrap_exp_βγ2)); + +/* Replace symbolic representations of β and γ expressions into their their +expansions */ +defrule(r_exp_β_val, β_expr(m_any1), β_expr_val(m_any1, expansion_order)); +defrule(r_exp_γ_val, γ_expr(m_any1), γ_expr_val(m_any1, expansion_order)); + + +/**/ + +freeof_many(expr, funcs) := apply(freeof, append(funcs, [expr])); + +/* Split product in factors containing any funcs and factors free of them */ +split_prod(expr, funcs) := block( + [deps_on_funcs : []], + /* If not a product, assume product with single factor */ + for f in funcs do block( + if atom(expr) then ( + if member(expr, funcs) then ( + push(expr, deps_on_funcs), + expr : 1), + return()) + elseif op(expr)="*" then block( + [split : partition(expr, f)], + expr : split[1], + push(split[2], deps_on_funcs)) + elseif not freeof_many(expr, funcs) then ( + push(expr, deps_on_funcs), + expr : 1, + return()) + else ( + return())), + [expr, apply("*", deps_on_funcs)]); + +/* Substitute functions by their series expansions and apply substitutions in +coeffs */ +__expand_term_subst(expr, coeffs) := block( + expr : apply1(expr, r_exp_β_val, r_exp_γ_val), + expr : subst(append(main_ansatz, coeffs), expr), + ev(expr, eval, diff, expand)); + +/* Substitute product of functions by result of mul_series */ +__expand_term_prod(expr, coeffs) := block( + [free, factors], + expr : wrap_exp(expr), + /* Split product into list of factors containing funcs and those independent + of funcs*/ + block([split : split_prod(expr, funcs)], + free : split[1], + if not atom(split[2]) and op(split[2])="*" then ( + factors : args(split[2])) + else ( + factors : [split[2]])), + /* Substitute funcs factors by expansion */ + factors : __expand_term_subst(factors, coeffs), + /* Total product using series optimization */ + free * mul_series(factors, expansion_order + lopow(free, r))); + +/* +Split product factors into list substitute series expansions, substitute +coefficients and finally multiply only upto expansion_order. +*/ +expand_term(expr, coeffs) := block( + [tmp], + /* TODO: Recursive call instead of error */ + if op(expr)="+" then ( + error("Term is itself a sum!")) + /* For single negation just call recursively */ + elseif op(expr)="-" then ( + if length(args(expr)) # 1 then ( + error("Unexpected '-' operator with more than 1 argument")), + expand(- expand_term(part(expr, 1), coeffs))) + elseif op(expr)="*" then ( + expand(__expand_term_prod(expr, coeffs))) + elseif op(expr)="/" then block( + [n, d, expansion_order:expansion_order - lopow(d, r)], + n : part(expr, 1), + d : part(expr, 2), + expand(expand_term(n, coeffs) / d)) + /* TODO: better handle exponents */ + else ( + apply1(__expand_term_subst(wrap_exp(expr), coeffs), r_simp_small))); + + +/* Solve main equations term by term */ + +/* TODO: Divide and conquer on term list for long expansions */ +/* coeffs: facts already known about series expansions */ + +/* Series expansion inside each instance of a value in the terms of an equation +*/ +expand_equation(expr, coeffs) := block( + [len, terms], + expr : rhs(expr) - lhs(expr), + 0 = if atom(expr) then ( + expr) + elseif op(expr)="+" then ( + apply("+", map(lambda([term], expand_term(term, coeffs)), args(expr)))) + else ( + expand_term(expr, coeffs)))$ + + +/* Misc */ + +/* For an equation of polynomials, returns the term by term list of equations of +the coefficients */ +coeff_equations(eq) := block([start, end], + eq : expand(eq), + start : lopow(eq, r), + end : hipow(eq, r), + reverse(makelist([r^i, coeff(eq, r, i)], i, start, end))); diff --git a/setup.mac b/setup.mac new file mode 100644 index 0000000..becfdc3 --- /dev/null +++ b/setup.mac @@ -0,0 +1,105 @@ +/* + * Setup + */ +linel: 120; + +derivabbrev: true; + +load("ctensor"); + +dim: 4; +cframe_flag: false; +ct_coords: [u, r, θ, ϕ]; + +/* Original Bondi U and V are replaced by A and B; V is reserved for potential */ +depends([A, B, β, γ], [u, r, θ]); + +/* +Bondi metric ansatz +NOTE: Signature is different than original Bondi paper +*/ +bondi_ansatz() := block([m], + m: zeromatrix(4, 4), + /* du du */ + m[1, 1] : -B * r^-1 * exp(2 * β) + A^2 * r^2 * exp(2 * γ), + /* du dr */ + m[1, 2] : -exp(2 * β), + m[2, 1] : -exp(2 * β), + /* du dθ */ + m[1, 3] : -A * r^2 * exp(2 * γ), + m[3, 1] : -A * r^2 * exp(2 * γ), + /* dθ dθ */ + m[3, 3] : r^2 * exp(2 * γ), + /* dϕ dϕ */ + m[4, 4] : r^2 * exp(- 2 * γ) * sin(θ)^2, + m)$ + +lg: bondi_ansatz()$ + +/* Load metric and compute inverse */ +cmetric()$ + +/* Force simplification of some expressions */ +ug[2,2] : expand(ug[2,2]),simp; + +/* mcs[i, j, k]: Christoffel symbols of the second kind (one upper index - k) */ +christof(false); + +/* ric[i, j]: R_ab Ricci tensor*/ +ricci(false); + +/* uric[i, j]: R^a_b Ricci tensor */ +uricci(false); + +/* Ricci scalar */ +ric_s : expand(scurvature()); + + +/* 2-sphere submanifold */ + +/* 2-sphere submanifold metric */ +h_AB[3,3] : exp(2*γ)$ +h_AB[4,4] : exp(-2*γ)*sin(θ)^2$ +h_AB[3,4] : 0$ +h_AB[4,3] : 0$ + +/* 2-sphere submanifold inverse metric */ +uh_AB[3,3] : exp(-2*γ)$ +uh_AB[4,4] : exp(2*γ)*sin(θ)^-2$ +uh_AB[3,4] : 0$ +uh_AB[4,3] : 0$ + +/* 2-sphere metric */ +q_AB[3,3] : 1$ +q_AB[4,4] : sin(θ)^2$ +q_AB[3,4] : 0$ +q_AB[4,3] : 0$ + +/* 2-sphere inverse metric */ +uq_AB[3,3] : 1$ +uq_AB[4,4] : sin(θ)^-2$ +uq_AB[3,4] : 0$ +uq_AB[4,3] : 0$ + +/* Christoffel symbols of 2-sphere */ +for i : 3 thru 4 do for j : 3 thru 4 do for k : 3 thru 4 do ( + h_mcs[i, j, k] : sum(uh_AB[k,C]*( + cdiff(h_AB[j, C], i) + + cdiff(h_AB[i, C], j) + + cdiff(h_AB[i, j], C))/2, C, 3, 4)); + +/* Riemann tensor of 2-sphere */ +for i : 3 thru 4 do for j : 3 thru 4 do for k : 3 thru 4 do for l : 3 thru 4 do ( + h_riem[i,j,k,l] : expand( + cdiff(h_mcs[i,k,l], j) + - cdiff(h_mcs[j,k,l], i) + + sum(h_mcs[i,k,n]*h_mcs[l,n,j] - h_mcs[j,k,n]*h_mcs[n,i,l], n, 3, 4))); + +/* Ricci tensor of 2-sphere */ +for i : 3 thru 4 do for j : 3 thru 4 do ( + h_ric[i,j] : expand(sum(h_riem[i,k,j,k], k, 3, 4))); + +/* Ricci scalar of 2-sphere */ +h_ric_s : expand(sum(sum(uh_AB[i,j]*h_ric[i,j], i, 3, 4), j, 3, 4)); + +print("Setup done!"); diff --git a/tensors.mac b/tensors.mac new file mode 100644 index 0000000..dea722a --- /dev/null +++ b/tensors.mac @@ -0,0 +1,21 @@ +/* Contract indices in expression */ + +/* NOTE: buildq works like splicing a list with quasiquote in lisp; escaped +items are defined at start of quote */ +contr(expr, var) ::= ( + if listp(var) and length(var) > 1 then + block([head, tail], + head : first(var), + tail : rest(var), + buildq([expr, head, tail], sum(contr(expr, tail), head, 1, 4))) + elseif listp(var) and length(var) = 1 then + block([head], + head : first(var), + buildq([expr, head], sum(expr, head, 1, 4))) + elseif listp(var) then + expr + else + buildq([expr, var], sum(expr, var, 1, 4)))$ + +/* Differentiate with respect to nth coordinate */ +cdiff(func, coord) := diff(func, ct_coords[coord])$ |