Differences between DSolve and NDSolve

Question

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}]

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}]

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

while the imaginary part of the solution gave this plot:

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?

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]

---

Addendum.

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).