my current project is numerically solve the Friedmann equation:

(a-dot / a)^2 = H_0^2 (sigma_r / a^4 + sigma_m / a^3 + sigma_lambda + simga_k / a^2)

What I did was split my equation into two differential equations. I wrote that

da/d(H_0t) = v

dv/d(H_0t) = simga_lambda*a - simga_r/a^3 - sigma_m / (2*a^2)

and then used the Runge-Kutta method to integrate this. You can check the math out and you will arrive at these two equations.

Here is my code to numerically solve the Friedmann equation:

(* Energy Densities *)
omegar =8.4*10^(-5);
omegam = 0.3;
omegalambda = 0.7;
omegacurve = 1 - omegar - omegam - omegalambda;

(* Time Interval *)
t0 = 10^(-6);
tf = 10;

(* Initial Conditions *)
a0 = Sqrt[ 2 * Sqrt[omegar] * t0]
y0 = Sqrt[ Sqrt[omegar] / (2*t0)]

eqn = {a'[t]==y[t], y'[t]==omegalambda*a[t] - omegar/ a[t]^3 - omegam /(2*a[t]^2), a[t0]== a0, y[t0]== y0};

sol = NDSolve[eqn, {a,y}, {t, t0, tf}]

LogLogPlot[a[t] /.First[sol], {t, t0, tf}, PlotRange->All]

Now, when I run this code, I receive the following error:

NDSolve::ndsz :  At  t  ==  0.000013790330904023948 , step size is effectively zero; singularity or stiff system suspected

Now, I can solve this problem when the step size is larger than 10^(-4), but not when it is smaller. I am using a different boundary condition for that case, but it shouldn't matter. I'm confused though on why this error keeps popping up and what I can do to solve this. What steps should I take for NDSolve?

  • 2
    $\begingroup$ Are you sure you've got your equations right? From your notation it looks like you're doing cosmology, but the equations you're using aren't equivalent to the Friedmann equation. Could you edit your question to include the ODE system you're trying to solve, along with a little context for it? (LaTeX code works on StackExchange, if you know how to use that.) $\endgroup$ Commented Jun 6, 2016 at 17:43
  • $\begingroup$ Yes, sure. What I did was that I took the time derivative of the Friedmann equation to avoid the square root if I didn't take it. I'll edit my post to account for it. $\endgroup$ Commented Jun 6, 2016 at 19:14

1 Answer 1


It looks like you're trying to integrate past the Big Crunch. Specifically, on a mathematical level, your equations are going to have trouble if $a$ goes to zero. You can use WhenEvent in your equations to detect when this is about to occur and stop the integration:

eqn = {a'[t] == y[t], 
       y'[t] == omegalambda*a[t] - omegar/a[t]^3 - omegam/(2*a[t]^2), 
       a[t0] == a0, y[t0] == y0, 
       WhenEvent[a[t] <= a0/10^5, tcrunch = t; "StopIntegration"]}

This tells Mathematica to stop integrating when $a(t)$ becomes less than $10^{-5} a(t_0)$. We then get the following result:

sol = NDSolve[eqn, {a, y}, {t, t0, tf}]
Plot[a[t] /. First[sol], {t, t0, Min[tcrunch, tf]}]

enter image description here

Note that the integration stops early, since the Universe has collapsed in on itself. (The precise time of the Big Crunch for this scenario is $t \approx 1.37903 \times 10^{-5}$.) The InterpolatingFunctions returned by NDSolve only extend from $t_0$ to $t_\text{crunch}.$

If you need to get closer to the Big Crunch, you can adjust the choice of 105 in the code above. When I brought it up to 106, I started getting the same errors about equation stiffness. (BTW, these types of errors are usually a sign that your solutions have blown up in some way.)

As to why this particular Universe is collapsing on itself, that's a question better suited for Physics StackExchange (appropriately recast, of course.)


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.