# Numerical solution for Friedmann equations

My problem today is to solve the Friedmann equations, for those who aren't familiar with them, here they are (in my specific case):

$$\left ( \frac{\dot{a}}{a^2} \right )^2 = \frac{\rho_1}{a^4} - \frac{\rho_2}{a^6}$$

So, my idea for solving this is to discretise in time, namely, write: $$\dot{a} = \pm \sqrt{\rho_1 - \frac{\rho_2}{a^2} }$$ And then $$\dot{a} = \frac{a_{i+1} - a_{i}}{dt}$$ Now, I want to start with $\dot{a} < 0$, and make it "bounce". Namely, if you solve these equations analytically, you get: $$a(t) = \sqrt{\rho_1 t^2 + \frac{\rho_2}{\rho_1}}$$ Meaning that $a(t)$ has a minimum value, at $t = 0$. However, using time discretisation one gets( in the phase when $\dot{a}<0$ $$a_{i+1} = a_i - dt \sqrt{\rho_1 - \frac{\rho_2}{a_i^2} }$$ Which at some point becomes complex.

This procedure has obviously some flaws, how can I correct it? I would like to write an algorithm that solves numerically using some sort of time discretisation, since I later will need to implement for a time variable $\rho_2(t)$. The algorithm must reproduce the analytical solution, with the feature that, once we have reached the minimum value for $a(t)$, it stops decreasing and starts increasing.

Any help is appreciated.

Rearranging your first equation you get $$\dot{a}^2 = \rho_1 - \frac{\rho_2}{a^2}$$ In this case you need only to give an initial condition $a(0)>\sqrt{\rho_2/\rho_1}$.
If you want to reach the bounce, integrate with a negative velocity forward in time (or vice versa), i.e. $\dot{a}=-\sqrt{\rho_1 -\rho_2/a^2}, \Delta t > 0$. Once you reach the bounce $a^2=\rho_2/\rho_1$ (or $\dot{a}=0$), switch the sign of the velocity so that your universe grows instead of shrinking. (Beware of NaNs from square-rooting a negative number!)
• @MrFermiMr Yes, sorry, I misread your first equation which is usually given in terms of $\dot{a}^2/a^2$ rather than $\dot{a}^2/a^4$ on the left hand side. Corrected.
• Yes, it's written in conformal time $dt = a(t') \; dt'$.. Actually i think that's precisely what happens: the solution goes and then, precisely at $a = \sqrt{\rho_1/\rho_2}$ the solution becomes complex. It is just a matter of telling mathematica to change sign when the complex part becomes non zero, and i should be good. Dec 25, 2014 at 10:13