I have a question about the difference in the solution between DSolve and NDSolve.

I want to solve the Friedmann equation of

\begin{align} 3\left(\frac{\dot{a}^2}{a^2}+\frac{k}{a^2}\right)=8\pi G\rho(a)\,, \end{align}

where $a(t)$ is the scale factor, $k=1,0,-1$, $G$ is Newton's gravitational constant, and $\rho$ is the energy density. For a Universe composed of radiation, I know that $\rho\sim a^{-4}$. For $k=1$ (and setting $G=1$), I can solve for $a(t)$ with Mathematica:

 DSolve[{3 (D[a[t], t]^2/a[t]^2 + 1/a[t]^2) == 8 Pi/a[t]^4, a[0] == 1}, a[t], t]

giving me a solution of:

{{a[t] -> Sqrt[3 - 2 Sqrt[3 (-3 + 8 \[Pi])] t - 3 t^2]/Sqrt[3]}, {a[t] -> Sqrt[3 + 2 Sqrt[3 (-3 + 8 \[Pi])] t - 3 t^2]/Sqrt[3]}}

Plotting the second solution, I get:

Plot[Sqrt[3 + 2 Sqrt[3 (-3 + 8 \[Pi])] t - 3 t^2]/Sqrt[2], {t, 0, 6}]

enter image description here

which is exactly what I want.

Now, solving the same equation with NDSolve:

sol = NDSolve[{3 (D[a[t], t]^2/a[t]^2 + 1/a[t]^2) == 8 Pi/a[t]^4, a[0] == 1}, a, {t, 0, 6}]

I get the following errors:

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

NDSolve::mxst: Maximum number of 799782 steps reached at the point t == 5.298557045191269.

Plotting the solution, I get

Plot[Evaluate[a[t] /. sol[[2]]], {t, 0, 5.3}]

enter image description here

We see that the first half of the plot looks similar to the analytical solution, but then nothing is plotted. I decided to plot the real part of the solution and I found

enter image description here

while the imaginary part of the solution gave this plot:

enter image description here

Why are there differences between DSolve and NDSolve? I tried adding the solution

 Method -> {"MethodOfLines", 
     "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 1}}

as suggested here, but it didn't fix the problem. I've also tried increasing MaxSteps to more than 100000 as well as changing the WorkingPrecision, PrecisionGoal, and AccuracyGoal. How do I make my numerical solution agree with my analytical solution?

  • $\begingroup$ But using DSolve, I get the error DSolve::bvnul: For some branches of the general solution, the given boundary conditions lead to an empty solution., don't you get it? I suspect you must inspect the validity of your analytic function. $\endgroup$ Nov 23, 2017 at 21:06

1 Answer 1


Both NDSolve and DSolve solve two equations:

Solve[{3 (D[a[t], t]^2/a[t]^2 + 1/a[t]^2) == 8 Pi/a[t]^4}, {a'[t]}]
{{Derivative[1][a][t] -> -(Sqrt[8 π - 3 a[t]^2]/(Sqrt[3] a[t]))},
 {Derivative[1][a][t] -> Sqrt[8 π - 3 a[t]^2]/(Sqrt[3] a[t])}}

The first message, NDSolve::ndsz is unimportant. It has to do with the solution to the first equation, which the OP ignores in both the DSolve and the NDSolve results.

The second comes up because a[t] approaches 8 π / 3 which causes the expression 8 π - 3 a[t]^2 under the radical to approximately vanish, perhaps even become negative due to numerical error at each discrete step.

A workaround is to differentiate the equation and choose for the additional initial condition, the value of the derivative a'[0] of the solution of interest.

Solve[{3 (D[a[t], t]^2/a[t]^2 + 1/a[t]^2) == 8 Pi/a[t]^4, a[t] == 1}, {a[t], a'[t]}]
{{a[t] -> 1, Derivative[1][a][t] -> -Sqrt[1/3 (-3 + 8 π)]},
 {a[t] -> 1, Derivative[1][a][t] -> Sqrt[1/3 (-3 + 8 π)]}}

The value for a'[0] is Sqrt[1/3 (-3 + 8 π)].

sol = NDSolve[{D[3 (D[a[t], t]^2/a[t]^2 + 1/a[t]^2) == 8 Pi/a[t]^4, t],
    a[0] == 1, a'[0] == Sqrt[1/3 (-3 + 8 π)]}, a, {t, 0, 6}];

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

ListLinePlot[a /. sol]

Mathematica graphics


I should probably have pointed out that a[t] == 8 π / 3 is a singular solution:

3 (D[a[t], t]^2/a[t]^2 + 1/a[t]^2) == 8 Pi/a[t]^4 /. 
  a[t] -> Sqrt[8 Pi/3] // Simplify
(*  Derivative[1][a][t] == 0  *)

When the NDSolve integration reaches a time t0 where a[t0] == 8 π / 3, there are two possible branches to follow. The one is shown in the DSolve solution (the same as mine above), and the other is the solution where a[t] == 8 π / 3 for all t > t0, which is approximately the NDSolve solution in the OP (looking at the real and imaginary parts).


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.