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$
5
  • $\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$
5
  • $\begingroup$ While this answer does provide the OP with a solution, I believe it should explain a bit more. $\endgroup$ – Alex Petrosyan Nov 29 '20 at 15:04
  • $\begingroup$ @AlexPetrosyan What do you need help with? It seems clear to me, but perhaps I am too familiar with NDSolve and Mathematica to see what I left out. $\endgroup$ – Michael E2 Nov 29 '20 at 15:41
  • $\begingroup$ Well, for starters, I think that the substitution semantics are a little obtuse. Why does your solution generate the substitution for both the function and its derivative, while OPs does not? What does the First @ do? $\endgroup$ – Alex Petrosyan Nov 29 '20 at 16:15
  • 1
    $\begingroup$ (1) It's explained in the answer: "which will replace a in both a[t] and a'[t] by the solution function." The solution does not return the derivative; it returns a function, namely an InterpolatingFunction. When replacing a in a'[t], the function is differentiated because that's what the prime means in Mathematica. The lesson here is that it is better to use NDSolve[..., a, {t, 0, 1}] than NDSolve[..., a[t], {t, 0, 1}].... $\endgroup$ – Michael E2 Nov 29 '20 at 16:43
  • 1
    $\begingroup$ (2) The @ is a way to denote function evaluation (on a single argument): f @ x is the same as f[x], only easier to write and to read, imo, because matching brackets is a painful chore to me. First[expr] is equivalent to expr[[1]] (or Part[expr, 1]) and returns the first element of the expression expr. Search for "First[" in the docs for NDSolve to see several uses, some of which are necessary. These are basic things that occur throughout Mathematica and don't need to be explained in every answer. $\endgroup$ – Michael E2 Nov 29 '20 at 16:45
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$
1
  • $\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.