Plotting the derivative of NDSolve divided by the solution

Question

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.

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$?

Answer 1

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

Answer 2

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}]

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