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?
DSolve
, I get the errorDSolve::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$