# Plotting the derivative of NDSolve divided by the solution

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

• 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? Jul 3 '19 at 1:13
• 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)$. Jul 3 '19 at 1:24
• Leave the equations alone. Just change the one place in front of {t, 0, 300*age} Jul 3 '19 at 1:31
• 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. Jul 3 '19 at 1:37
• I finally understand, could've stared at the code for hours and would never have noticed such a tiny detail. Thank you very much. Jul 3 '19 at 1:41

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 == 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 == 10^-10,
b'[t] = 1/a[t], b == 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 == 10^-10,
b'[t] == 1/a[t], b == 0}, {a, b, dloga}, {t, 0, 300*age}]

ListLogLogPlot[# /. scale, PlotLabel -> #, Joined -> True] & /@ {a, b, dloga} //
GraphicsRow • While this answer does provide the OP with a solution, I believe it should explain a bit more. Nov 29 '20 at 15:04
• @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. Nov 29 '20 at 15:41
• 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? Nov 29 '20 at 16:15
• (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}].... Nov 29 '20 at 16:43
• (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[] (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. Nov 29 '20 at 16:45

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 == 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 == 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 == 10^-10}, a, {t, 0, 300*age}]
Plot[Evaluate[a'[t]/a[t] /. scaleb[[1, 1]]], {t, 0,  age}]


same picture

• 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)}$? Jul 3 '19 at 1:38