4
$\begingroup$

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]

enter image description here

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

3 Answers 3

5
$\begingroup$

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

enter image description here

$\endgroup$
3
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – Rudyard
    Aug 23, 2022 at 14:32
1
$\begingroup$

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]

enter image description here

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.