Pertubation theory in sagemath

Question

I have seen the documentation on how to truncate polynomials using sage but I am stuck as to how I can actually apply this in my work...

I am currently trying to find the...say Ricci tensor for a perturbed metric tensor. Adding the small amount of free scalars (in this case, 2) increases the computation time by a lot because Sage is currently attempting to compute the solution EXACTLY. The problem is, I don't NEED the exact solution, only a 'first order' solution. When doing these calculations by hand, it is easy to throw away terms you don't need because they are quadratic order or higher, but I want to be able to do this using sage... Basically, I want to know if I can tell sage to throw away terms of higher order during the calculation so the computation time decreases.

PS: I have attempted to use symbolic rings to use $\mathcal{O}$ notation and this works perfectly on functions but this does not work when attempting to merge this with a Lorentzian metric.

For example, if my perturbations are $\psi$ and $\zeta$ I don't want sage to consider terms that are overall $\zeta$$\psi$ or $\zeta^2$ or $\psi^2$ etc.

Answer 1

You can accomplish this by introducing a small parameter $\epsilon$ and using set_calc_order. Here's an example of computing the Ricci scalar for a perturbed Schwartzchild metric to $\mathcal{O}(\epsilon^2)$.

First set up the manifold in sage:

%display latex
M = Manifold(4, 'M', latex_name=r'\mathcal{M}', structure='Lorentzian')
chart.<t,r,th,ph> = M.chart(r"t r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\phi")
m = var('m')
assume(m > 0)

Introduce a symbolic infinitesimal parameter:

eps = var('epsilon')

Define the metric (or scalar field or whatever else you want) using eps:

g = M.metric()
g[0,0] = -(1-2*m/r) + eps
g[1,1] = 1/(1-2*m/r) + eps
g[2,2] = r^2
g[3,3] = (r*sin(th))^2
g.display()

$$\newcommand{\Bold}[1]{\mathbf{#1}}g = \left( \epsilon + \frac{2 \, m}{r} - 1 \right) \mathrm{d} t\otimes \mathrm{d} t + \left( \epsilon - \frac{1}{\frac{2 \, m}{r} - 1} \right) \mathrm{d} r\otimes \mathrm{d} r + r^{2} \mathrm{d} {\theta}\otimes \mathrm{d} {\theta} + r^{2} \sin\left({\theta}\right)^{2} \mathrm{d} {\phi}\otimes \mathrm{d} {\phi}$$

The full Ricci scalar without truncation is:

ricci = g.ricci_scalar()
ricci.expr().factor()

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{2 \, {\left(4 \, \epsilon^{3} m^{2} r^{2} - 4 \, \epsilon^{3} m r^{3} + \epsilon^{3} r^{4} + 16 \, \epsilon^{2} m^{3} r - 28 \, \epsilon^{2} m^{2} r^{2} + 12 \, \epsilon^{2} m r^{3} - \epsilon^{2} r^{4} + 16 \, \epsilon m^{4} - 48 \, \epsilon m^{3} r + 32 \, \epsilon m^{2} r^{2} - 4 \, \epsilon m r^{3} - \epsilon r^{4} - 12 \, m^{4} + 12 \, m^{3} r + 2 \, m^{2} r^{2} - 4 \, m r^{3} + r^{4}\right)} \epsilon}{{\left(2 \, \epsilon m - \epsilon r - r\right)}^{2} {\left(\epsilon r + 2 \, m - r\right)}^{2} r^{2}}$$

Set the order of $\epsilon$ with set_calc_order. Now computations are truncated to $\mathcal O(\epsilon^2)$ automatically:

g.set_calc_order(eps, 2)

ricci = g.ricci_scalar()
ricci.expr().factor()

$$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{2 \, {\left(96 \, \epsilon m^{6} - 224 \, \epsilon m^{5} r + 192 \, \epsilon m^{4} r^{2} - 64 \, \epsilon m^{3} r^{3} + 6 \, \epsilon m r^{5} - \epsilon r^{6} + 24 \, m^{5} r - 36 \, m^{4} r^{2} + 8 \, m^{3} r^{3} + 10 \, m^{2} r^{4} - 6 \, m r^{5} + r^{6}\right)} \epsilon}{{\left(2 \, m - r\right)}^{3} r^{5}}$$