Solve the perturbed Boltzmann Einstein equations in cosmology using RK4 method

Question

I'm trying to learn to numerically solve the perturbed Boltzmann-Einstein equations in cosmology using the RK4 method. These are the equations:

\begin{align} \dot{\Theta}_{r,0}+k\Theta_{r,1}&=-\dot{\Phi} \\ \dot{\Theta}_{r,1}+\frac{k}{3}\Theta_{r,0} & =\frac{-k}{3}\Phi \\ \dot{\delta}+ikv &= -3\dot{\Phi} \\ \dot{v}+\frac{\dot{a}}{a}v &= ik\Phi \\ \dot{\Phi}&=\frac{1}{3\dot{a}}\frac{3H_{0}^{2}}{2}\left(\Omega_{m}\delta+4\Omega_{r}\Theta_{r,0}a^{-1}\right)-ak^{2}\Phi-\frac{\dot{a}}{a}\Phi \end{align}

How do I go about this? I am feeling lost to even start.

Answer 1

The only difference between your system of ODEs to the single variable ODE is that rather than returning a single value for each step, you are returning a list/vector/array (depending on your programming language) with 1 element for each independent parameter.

As an example, consider this arbitrary system of two coupled ODEs, $$\frac{\mathrm{d}x}{\mathrm{d}t} = f(t,x,y);\quad\frac{\mathrm{d}y}{\mathrm{d}t}=g(t,x,y)\Longrightarrow \frac{\mathrm{d}}{\mathrm{d}t}\left[\begin{array}{c}x \\ y\end{array}\right]=\left[\begin{array}{c}f(t,x,y) \\ g(t,x,y)\end{array}\right]$$ with $x(0)=x_0$ and $y(0)=y_0$. This is solved mostly the same as the one-parameter case (assuming 0-indexed language, some languages out there are 1-indexed):

k1 = [h * f(t, x, y), h * g(t, x, y)]
k2 = [h * f(t+h/2, x + k1[0]/2, y + k1[1]/2), h * g(t+h/2, x + k1[0]/2, y + k1[1]/2)]
k3 = [h * f(t+h/2, x + k2[0]/2, y + k2[1]/2), h * g(t+h/2, x + k2[0]/2, y + k2[1]/2)]
k4 = [h * f(t+h, x + k3[0], y + k3[1]), h * g(t+h, x + k3[0], y + k3[1])]
x = x + (k1[0] + 2 * (k2[0] + k3[0]) + k4[0]) / 6
y = y + (k1[1] + 2 * (k2[1] + k3[1]) + k4[1]) / 6
t = t + h

In your case, you just need to write the system as, $$\frac{\mathrm{d}}{\mathrm{d}t}\left[\begin{array}{c}\Theta_{r,0} \\ \Theta_{r,1} \\ \delta \\ v \\ \Phi\end{array}\right]=\left[\begin{array}{c}f_1(t,\,\Theta_{r,0},\,\Theta_{r,1},\,\delta,\,v,\,\Phi,\,\dot{\Phi}) \\ f_2(t,\,\Theta_{r,0},\,\Theta_{r,1},\,\delta,\,v,\,\Phi,\,\dot{\Phi}) \\ f_3(t,\,\Theta_{r,0},\,\Theta_{r,1},\,\delta,\,v,\,\Phi,\,\dot{\Phi}) \\ f_4(t,\,\Theta_{r,0},\,\Theta_{r,1},\,\delta,\,v,\,\Phi,\,\dot{\Phi}) \\ f_5(t,\,\Theta_{r,0},\,\Theta_{r,1},\,\delta,\,v,\,\Phi,\,\dot{\Phi})\end{array}\right]$$ (and if you need to add $\dot{a}=f_6(\cdots)$, that should be straight-forward also), then it can be solved much the same as the 2-variable case outlined above.