substituting oscillatory solution with an averaged version efficiently

Question

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]

Answer 1

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

Answer 2

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.

Answer 3

My second answer. Numerical solution.

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