3
$\begingroup$

I'm trying to solve the Friedmann equation as is given by:

$$ H^{2}=\Big(\frac{\dot{a}}{a}\Big)^{2}=H_{0}^{2}\big(\Omega_{m,0}a^{-3}+\Omega_{r,0}a^{-4}+\Omega_{\Lambda,0}a^{2.3}+\Omega_{k}a^{-2}\big),$$

with the following values of certain parameters $$\Omega_{m,0}=0.3, \,\Omega_{r,0}=0.001, \,\Omega_{\Lambda,0}=0.6, \, \Omega_{k}=1-\Omega_{r,0}-\Omega_{m,0}-\Omega_{\Lambda, 0} $$

and $w=-0.9$ where $w$ determines the peculiar $2.3$ from the equation above. To solve, naturally, I used Mathematica's NDSolve function:

m = 0.3; 
r = 0.001;
d = 0.6;
k = 1 - m - r - d;
Hubble = 2.33*10^-18; age = 1/Hubble; 

scale = NDSolve[{a'[t]/a[t] == 
    Hubble*Sqrt[m*a[t]^-3 + r*a[t]^-4 + d*a[t]^-0.6 + k*a[t]^-2], 
   a[0] == 10^-10}, a[t], {t, 0, 300*age}]

Plotting this with Plot[Evaluate[{a[t], a'[t]} /. scale2], {t, 0, age}] I get the bottom left picture. Trying to draw the derivative of this function results in failure (following the GetHelp section and the numerous answers to variants of this question on this board did not help) which is shown on the bottom right. However, my endgoal isn't just to draw the derivative, I want to plot how the Hubble parameter $H(t)=\frac{\dot{a}}{a}$ varies with time.

enter image description here

My attempts at solving:

Drawing just the derivative: 

Plot[Evaluate[{a'[t]} /. scale2], {t, 0, age}, PlotStyle -> Blue]

Trying to plot the Hubble parameter:

Plot[Evaluate[{a'[t]/a[t]} /. scale2], {t, 0, age}]

which results in the same "empty" picture when trying to plot the derivative by itself. The blue above was added so I might hopefully notice it somewhere on the graph but to no avail. I've also tried using and ignoring Evaluate but no option worked. Also, is it possible to work with $a(t)$ now inside an integral, i.e. integrate an equation of the form $\int f\big(a(t)\big)\,dt$?

| improve this question | |
$\endgroup$
  • $\begingroup$ First solve for the function a not the value a[t] and then remove the extra {} with First: scale = First@NDSolve[..., a, {t, 0, 300*age}]. Does it work? $\endgroup$ – Michael E2 Jul 3 '19 at 1:13
  • $\begingroup$ I'm not quite sure what you mean. Do I swap a[t] with a everywhere it shows up in scale or? If I've understood correctly, I've tried that but then nothing plots at all, not even the solution $a(t)$. $\endgroup$ – Kandrax Jul 3 '19 at 1:24
  • $\begingroup$ Leave the equations alone. Just change the one place in front of {t, 0, 300*age} $\endgroup$ – Michael E2 Jul 3 '19 at 1:31
  • 1
    $\begingroup$ The solution in your code is of the form {{a[t] -> ...}}, which replaces only a[t], not a'[t]. If you make the two changes, it will become {a -> ...}, which will replace a in both a[t] and a'[t] by the solution function. $\endgroup$ – Michael E2 Jul 3 '19 at 1:37
  • $\begingroup$ I finally understand, could've stared at the code for hours and would never have noticed such a tiny detail. Thank you very much. $\endgroup$ – Kandrax Jul 3 '19 at 1:41
3
$\begingroup$

This should work:

scale = First@NDSolve[{a'[t]/a[t] == 
    Hubble*Sqrt[m*a[t]^-3 + r*a[t]^-4 + d*a[t]^-0.6 + k*a[t]^-2], 
   a[0] == 10^-10}, a, {t, 0, 300*age}]

Plot[Evaluate[{a'[t]/a[t]} /. scale], {t, 0, age}]

The solution in the OP's code is of the form {{a[t] -> ...}}, which replaces only a[t], not a'[t]. With the changes above, the solution scale becomes {a -> ...}, which will replace a in both a[t] and a'[t] by the solution function.

To get the (indefinite) integral of 1/a[t], add it to the system:

scale = First@NDSolve[{a'[t]/a[t] == 
    Hubble*Sqrt[m*a[t]^-3 + r*a[t]^-4 + d*a[t]^-0.6 + k*a[t]^-2], 
   a[0] == 10^-10,
   b'[t] = 1/a[t], b[0] == 0}, {a, b}, {t, 0, 300*age}]

Then the integral may be obtained with b /. scale.

Update

You can get a'[t]/a[t] (= D[Log[a[t]],t]) right from the system, just like the integral:

scale = First@ NDSolve[{
    dloga[t] == a'[t]/a[t] == 
     Hubble*Sqrt[m*a[t]^-3 + r*a[t]^-4 + d*a[t]^-0.6 + k*a[t]^-2],
    a[0] == 10^-10,
    b'[t] == 1/a[t], b[0] == 0}, {a, b, dloga}, {t, 0, 300*age}]

ListLogLogPlot[# /. scale, PlotLabel -> #, Joined -> True] & /@ {a, b, dloga} //
  GraphicsRow

enter image description here

| improve this answer | |
$\endgroup$
2
$\begingroup$

You can solve for a'[t]/a[t] using NDSolveValue:

aprimeovera =  NDSolveValue[{a'[t]/a[t] == Hubble*Sqrt[m*a[t]^-3 + 
   r*a[t]^-4 + d*a[t]^-0.6 + k*a[t]^-2],  a[0] == 10^-10}, 
  a'[t]/a[t], {t, 0, 300*age}]

Plot[aprimeovera, {t, 0, age}]

enter image description here

Update:

scale = NDSolve[{a'[t]/a[t] == 
    Hubble*Sqrt[m*a[t]^-3 + r*a[t]^-4 + d*a[t]^-0.6 + k*a[t]^-2], 
   a[0] == 10^-10}, a[t], {t, 0, 300*age}]

Convert solution to a pure function:

fa = Function[t, Evaluate[a[t] /. scale[[1, 1]]]];
Plot[Evaluate[fa'[t]/fa[t]], {t, 0,  age}]

same picture

Or, as MichaelE2 suggested, use a in the second argument of NDSolve to get a pure function:

scaleb = NDSolve[{a'[t]/a[t] == 
    Hubble*Sqrt[m*a[t]^-3 + r*a[t]^-4 + d*a[t]^-0.6 + k*a[t]^-2], 
   a[0] == 10^-10}, a, {t, 0, 300*age}]
Plot[Evaluate[a'[t]/a[t] /. scaleb[[1, 1]]], {t, 0,  age}]

same picture

| improve this answer | |
$\endgroup$
  • $\begingroup$ Simple yet effective, thank you. Any ideas as to why Mathematica refuses to do it in the "standard" way? Also, any suggestions on how I might go about solving integrals that include $a(t)$, say, $\int \frac{1}{a(t)}$? $\endgroup$ – Kandrax Jul 3 '19 at 1:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.