Ask Your Question

Can't substitute H=dot(a)/a in SageMath

asked 2020-02-14 21:55:09 +0200

cduston gravatar image

Hello all, I'm working with the Friedman equations, and I've gotten them down to the common form presented in terms of the scale factor. MWE incoming:

M = Manifold(4, 'M', structure='Lorentzian')
fr.<t,r,th,ph> = M.chart(r't r:[0,+oo) th:[0,pi]:\theta ph:[0,2*pi):\phi')
var('G, Lambda, k', domain='real')
a = M.scalar_field(function('a')(t), name='a')
rho = M.scalar_field(function('rho')(t), name='rho')
p = M.scalar_field(function('p')(t), name='p')
g = M.metric()
g[0,0] = -1
g[1,1] = a*a/(1 - k*r^2)
g[2,2] = a*a*r^2
g[3,3] = a*a*(r*sin(th))^2
nabla = g.connection()
Ricci = nabla.ricci()
Ricci_scalar = g.ricci_scalar()
E1 = Ricci - Ricci_scalar/2*g + Lambda*g
print("First Friedmann equation:\n")
E1[0,0].expr().expand() == 0

With the result being

-Lambda + 3*diff(a(t), t)^2/a(t)^2 + 3*k/a(t)^2 == 0

Which is great, but I really want this to be in terms of the Hubble paramter $H=\dot{a}/a$. I can't make that substitution happen:


just spits out the same thing,

-Lambda + 3*diff(a(t), t)^2/a(t)^2 + 3*k/a(t)^2 == 0

Any ideas?

edit retag flag offensive close merge delete

2 Answers

Sort by ยป oldest newest most voted

answered 2020-02-14 22:21:39 +0200

eric_g gravatar image

Assuming that eq1 is E1[0,0].expr().expand() == 0, what about

sage: eq1.subs({diff(a.expr(),t): h*a.expr()})
3*h^2 - Lambda + 3*k/a(t)^2 == 0

Alternatively, you may want to keep h being a function of t:

sage: eq1.subs({diff(a.expr(),t): function('h')(t)*a.expr()})
3*h(t)^2 - Lambda + 3*k/a(t)^2 == 0
edit flag offensive delete link more

answered 2020-02-14 23:26:28 +0200

rburing gravatar image

I hope you appreciate that this is not such a simple operation to perform in general.

The simplest way is probably to replace $\dot{a}$ by $Ha$.

We can use the class I gave as an answer to Substituting a differential equation into an expression:

sage: eq1 = E1[0,0].expr().expand() == 0
sage: var('h')
sage: pde = diff(a.expr(), t) == h*a.expr()
sage: s = SubstituteEvolutionaryPDE(pde, t)
sage: s(eq1)
3*h^2 - Lambda + 3*k/a(t)^2 == 0

This also works when there are second or mixed derivatives, e.g.:

sage: eq2 = diff(a.expr(), t, t) + 13 == 0
sage: s(eq2)
h^2*a(t) + 13 == 0
edit flag offensive delete link more


Well, since this is my major frustration with Sage, I know that it is not a simple operation, but I can't for the life of me understand why it should be. Although your answer does the job, I'm going to select the other one as it doesn't require a new class definition. Although please, disagree and convince me!

cduston gravatar imagecduston ( 2020-02-19 23:01:01 +0200 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower


Asked: 2020-02-14 21:55:09 +0200

Seen: 80 times

Last updated: Feb 14 '20