# substituting oscillatory solution with an averaged version efficiently

I am numerically solving the below ODE which has a highly oscillatory solution. This is fine.

Mpl = 10^19; \[CapitalLambda]qcd = 2/10;
m = 1;
ti = 2 10^12;
tf = 4 10^14;
eqns={a''[t] + 3/(2 t) a'[t] + m^2 ((2 \[CapitalLambda]qcd^2 t)/Mpl )^4 a[t] == 0, a[ti] == 1, a'[ti] == 0};
st = NDSolve[eqns, a[t], {t, ti, tf}, WorkingPrecision -> 10];
En[t_] := (t^(7/2) a[t]^2)/(3.48 10^46) /. st
LogLinearPlot[En[t], {t, ti, tf}, PlotRange -> All]


However I would like to substitute what I'm plotting (the function En[t] in code) with an 'averaged' version. I've gone about this the explcit way - see below - by defining a new function which at every point averages over many cycles but plotting this takes forever. There must be a more efficient way? If not in defining the new function, at least in plotting it...

s\[Tau] = st /. t -> \[Tau];
EnAvg[t_] := NIntegrate[(a[\[Tau]] /. s\[Tau])^2/(0.2 t), {\[Tau], 0.9 t , 1.1 t}][[1]]
(* Takes forever *)
LogLinearPlot[EnAvg[t], {t, ti, tf}, PlotRange -> All]

• Side question: Why WorkingPrecision -> 10? It's slower and less precise that using machine precision. Is that what you want? Aug 23, 2022 at 12:50
• You can use the NDSolve method in this answer to integrate a[t]^2 (instead of the inverse function), using NDSolveIterate[state, #] & /@ Flatten[a[t] /. st /. t -> "Grid"] for the Iterate stage. Or you can simultaneously integrate a[t]^2 in the original NDSolve call. Aug 23, 2022 at 13:19

Here's one way, among many...

Grab the points generated for your plot (I find this trick repeatedly useful in working with Mathematica)

plot = LogLinearPlot[En[t], {t, ti, tf}, PlotRange -> All];
points = Cases[plot, Line[data_] :> data, Infinity][[1]];

points[[1;;3]]
(* {{28.3242, 0.000325107},
{28.3258, 0.000326961},
{28.3274, 0.000328826}}
*)


Since the Plot function is adaptive, the X-axis values aren't sampled on a fixed grid, so use MovingMap to smooth the data. You get to pick the window and the smoothing function. Two obvious choices are Mean and Median. The third input is the window size. The None will ensure the same number of points in the smoothed data as the original, but the smoothing window "overhangs" the data at the end.

smooth1 = MovingMap[Mean, points, .3, None];
smooth2 = MovingMap[Mean, points, 3, None];
smooth3 = MovingMap[Median, points, .3, None];
smooth4 = MovingMap[Median, points, 3, None];

ListLinePlot[{smooth1, smooth2, smooth3, smooth4}]


You can solve this diffequation anylytically with DSolve.

Mpl = 10^19; \[CapitalLambda]qcd = 2/10;
m = 1;
ti = 2 10^12;
tf = 4 10^14;
eqns = {a''[t] + 3/(2 t) a'[t] +
m^2 ((2 \[CapitalLambda]qcd^2 t)/Mpl)^4 a[t] == 0, a[ti] == 1,
a'[ti] == 0} // Rationalize[#, 0] &;

eqns2 = Subtract @@@ eqns // Together // Numerator;

asol = a /. First@DSolve[Thread[eqns2 == 0], a, t]

En[t_] = (t^(7/2) asol[t]^2)/(348 10^44) // FullSimplify

(*   ((2 + Sqrt[3]) (t^3)^(
5/6) (\[Pi] (t^3)^(1/6)
BesselJ[-(13/12), 8/46875] BesselJ[1/12, t^3/
46875000000000000000000000000000000000000] + \[Pi] Sqrt[t]
BesselJ[-(1/12), t^3/
46875000000000000000000000000000000000000] BesselJ[13/12, 8/
46875])^2)/(597381591796875000000000000000000000000000000000 \
Sqrt[2] Sqrt[t])   *)

Plot[En[t], {t, ti, tf }, WorkingPrecision -> 50, PlotPoints -> 100,
PlotRange -> All]


Plot is different from your plot, since your plot suffers of workingprecision for higher t. Now regard En2[t] to see the two only osccilarory BesselJ terms, which oscillate to zero for high t.

rules = {BesselJ[1/12, t^3/
46875000000000000000000000000000000000000] -> bes1,
BesselJ[-(1/12), t^3/46875000000000000000000000000000000000000] ->
bes2};

En2[t_] = En[t] /. rules

(*   ((2 + Sqrt[3]) (t^3)^(
5/6) (bes1 \[Pi] (t^3)^(1/6) BesselJ[-(13/12), 8/46875] +
bes2 \[Pi] Sqrt[t]
BesselJ[13/12, 8/
46875])^2)/(597381591796875000000000000000000000000000000000 \
Sqrt[2] Sqrt[t])   *)

LogLinearPlot[
BesselJ[1/12, t^3/46875000000000000000000000000000000000000], {t, ti,
tf}, PlotRange -> All, GridLines -> Automatic,
WorkingPrecision -> 50, PlotPoints -> 100]


I leave it to you to extract behavior for higher t, averaging oscillaroty terms.

• Thank you. I know in this particular simple case one can solve analytically but I am actually interested in more complicated versions where there can only be a numerical handle. Aug 23, 2022 at 14:32

Use NDSolve to integrate En in the same run as calculating a.

Mpl = 10^19; \[CapitalLambda]qcd = 2/10;
m = 1;
ti = 2 10^12;
tf = 4 10^14;
eqns = {a''[t] + 3/(2 t) a'[t] +
m^2 ((2 \[CapitalLambda]qcd^2 t)/Mpl)^4 a[t] == 0, a[ti] == 1,
a'[ti] == 0} // Rationalize[#, 0] &;

eqns2 = Subtract @@@ eqns // Together // Numerator;

nsol = First@
NDSolve[Join[
Thread[eqns2 == 0], {Enint'[t] == (t^(7/2) a[t]^2)/(348 10^44),
Enint[ti] == 0}], {a, Enint}, {t, ti, tf}, MaxSteps -> 5 10^5,
WorkingPrecision -> 80, AccuracyGoal -> 10]

Plot[a[t] /. nsol, {t, ti, tf}, WorkingPrecision -> 70,
PlotRange -> All, PlotPoints -> 100]

Plot[Enint[t] /. nsol, {t, ti, tf}, WorkingPrecision -> 70,
PlotRange -> All, PlotPoints -> 100]


You see, oscillations do not compensate to constant value, but continuous integral of En is growing. But please check, whether accuracy is high enough.

Get En[t] by simply differetiating the integral Enint[t].

Plot[Enint'[t] /. nsol, {t, ti, tf}, WorkingPrecision -> 70,
PlotRange -> All, PlotPoints -> 100]


You see, En[t] has a mean value of about 2.5, which is the gradient/slope in the picture of Enint[t].

(Enint[4 10^14] - Enint[2 10^14])/(2 10^14) /. nsol
`

(* 2.653850693851505840381830 *)

More exact: Linear fit