Funky behaviour of derivative
The problem is that the default interpolating functions returned by NDSolve
are cubic Hermite interpolations. The first derivative will be continuous but (most likely) with singularities at the steps (values of t
for each step). The second derivative will be pretty much useless. The fix is to use InterpolationOrder -> All
(and hope for the best).
With the InterpolationOrder -> All
, the interpolation with the default LSODA method is primarily by local series, mainly of order 4 or 5. It seems impossible to completely avoid cubic Hermite interpolation with the LSODA method, at least with this ODE. However, if we use Method -> "ExplicitRungeKutta"
and set up the ODE system in a nicer way, we can get an interpolation completely in terms of Chebyshev series order 8 or 9, which are smooth enough for plotting the OP's functions.
One significant change is differentiating the first ODE and explicitly solving for H'[t]
and x''[t]
for NDSolve
. For some reason NDSolve[]
had trouble solving for the derivatives. Also, as @bbgodfrey pointed out in a comment, a[t]
does not have a initial value; in the OP's code, NDSolve[]
arbitrarily assign an initial value of 1
.
(* use OP's parameters and other definitions *)
ics = {xp - Sqrt[2/3]*mp*Log[ti] == x[ti], -Sqrt[2/3]*mp/ti == x'[ti],
1/(3*ti) == H[ti],
(*H'[ti] == 1/3,*) (* problematic & seems to be wrong -- omit *)
a[ti] == 1}; (* this could be anything -- OP should specify *)
ode = {Equal @@@ First@Solve[{
D[H[t]^2 - 8/3 π (x[t]^2/(80000000000 π) + 1/2 x'[t]^2) == 0, t],
x[t]/(40000000000 π) + 3 H[t] x'[t] + x''[t] == 0},
{H'[t], x''[t]}],
a'[t]/a[t] == H[t]};
solution = NDSolve[{ode, ics},
{H, x, a}, {t, 0.1, 10^2},
InterpolationOrder -> All, Method -> "ExplicitRungeKutta",
AccuracyGoal -> Infinity, MaxSteps -> Infinity];
z[t_] = a[t]*x'[t]/H[t] /. solution;
w[t_] = (D[z[t], {t, 2}]*(a[t]^2 /. solution) +
z'[t]*(a'[t] /. solution)*(a[t] /. solution))/z[t];
Then the plot looks like this:
Plot[{z[t], z'[t], z''[t], w[t]}, {t, 1, 10}, PlotLegends -> "Expressions"]
Note: You can check the principal interpolation method with
x["InterpolationMethod"] /. solution
(* {"Chebyshev"} *)
It may be of interest to present a solution requiring minimal changes to the code in the question, namely adding a[ti] == 1
(so that the code runs on Ver 11.0.1; DSolve
chooses a[ti] == 1
by default in earlier versions) and InterpolationOrder -> 5
(to eliminate the spiky behavior of some of the curves in the question, as I noted in an early comment). With these changes only, the curves are the same as in the solution by Michael E2, apart from a small spike in w
at about t == 3.4
(which can be eliminated with the Plot
option PlotPoints -> 20
. Curves of the three functions obtained directly from NDSolve
are,
LogLogPlot[(H /. Flatten@solution)[t], {t, 1, 10^8}]
LogLogPlot[(a /. Flatten@solution)[t], {t, 1, 10^8}]
Plot[(x /. Flatten@solution)[t], {t, 1, 10^8}, PlotRange -> {-.05, .05}]
The small PlotRange
is used to make visible the late-time oscillations.