3
$\begingroup$

I am trying to solve a system of second order differential equation, over time, which has periods of fast changing profile, and I need to evaluate the equation around that transition. I have to run it over a scan of points, thus time per point needs to be below O(1) second. The trouble is that for some of the values, the ODE doesn't wanna solve and spits "Indeterminate" or "Complex Infinity". I tried using MaxStepSize, but that only leads to bogus values as results. I also tried to increase the time interval over which my equation is being evaluated, and although it is giving sensible answers, it takes O(1000) seconds for a single point, making it untenable.

I am attaching the code for one specific point, any suggestions would be greatly appreciated. Let me know if there are any further clarifications needed.


(*VALUES FOR VARIABLES USED*)

fa = 1.13*10^13;
Rf = 0.0104;
\[Theta] = 1;  
\[Phi] = 1;
Mpl = 12*10^18; 
\[CapitalLambda]qcd = 2/10;
f\[CurlyPhi] = Rf* fa;
ma$t[t_] = 
 If[t < Mpl/(2 \[CapitalLambda]qcd^2), ma (\[CapitalLambda]qcd/Sqrt[Mpl/(2 t)])^4,  ma]; 

t1 = (Mpl^2/(8 ma \[CapitalLambda]qcd^4))^(1/3);
ma =  57*10^-16 (10^12/fa); 

Rm = Rf^4 (3/10 87/10 24/10 1/1000 (fa/10^12)^(3/2))/(11/100);
m\[CurlyPhi] = Rm ma;

(* Mass matrix eigenvalues *)

mH$t[t_] = Sqrt[1/2 ma$t[t]^2 (1 + Rm^2 ma^2/ma$t[t]^2 (1 + Rf^2 + 
        Sqrt[(1 - Rf^2 - 1/Rm^2 ma$t[t]^2/ma^2)^2 + 4 Rf^2]))];

mL$t[t_] = Sqrt[1/2 ma$t[t]^2 (1 + Rm^2 ma^2/ma$t[t]^2 (1 + Rf^2 - 
        Sqrt[(1 - Rf^2 - 1/Rm^2 ma$t[t]^2/ma^2)^2 + 4 Rf^2]))];

(*Initial and final solving time, increasing the final times makes the equation take forever*)
ti = t1/1000; 
tf = t1 250;

(*DIFFERENTIAL EQUATION*)

eqns = {a''[t] + 3/(2 t) a'[t] + ma$t[t]^2 a[t] + m\[CurlyPhi]^2 Rf^2 a[t] + 
       m[CurlyPhi]^2 Rf\[CurlyPhi][t] == 0, \[CurlyPhi]''[t] + 3/(2 t) \[CurlyPhi]'[t] + 
       m\[CurlyPhi]^2 \[CurlyPhi][t] + m\[CurlyPhi]^2 Rf a[t] == 0, a'[ti] == 0,           
       \[CurlyPhi]'[ti] == 0, a[ti] == \[Theta] fa, \[CurlyPhi][ti] ==                   
       \[Phi]f[CurlyPhi]};

{{sol1, sol2}} = NDSolve[eqns, {a[t], \[CurlyPhi][t]}, {t, ti, tf}, MaxSteps -> Infinity, 
                 PrecisionGoal -> 6];

UPDATE : Thanks a lot to Alex, and following his solution involving dimensionless quantities and NDSolve method tricks, it speeded up the evaluation, and also realizing that I had to solve until t0 at max, I used "ADAMS" method. But for some points, I realized that using "BDF" actually resolved the issue better, so would like to understand why's that.

The query can be thought in a generic way : which Methods and tactics are best used for solving 2nd order ODEs, which have solutions which start at some initial value, slowly roll towards it's minima and then oscillate around it.

As an example, am putting in the updated code:


\[Theta] = 1; 
\[Phi] = 1;

fa = 2.1544346900318835855901063672841`10.*^12;

ma = 57 10^-16 (10^12/fa);  (* GeV axion mass *)
t1 = (Mpl^2/(
  8 ma \[CapitalLambda]qcd^4))^(
 1/3); (* GeV^-1 time when ma(T) = H(T) *)

t0 = Mpl/(2 \[CapitalLambda]qcd^2)/t1;

Rf = 0.6256495158757365`

fphi = f\[CurlyPhi] = Rf* fa;

Rm =0.0116328;

mphi = m\[CurlyPhi] = Rm ma;


m0[t_] := ma (\[CapitalLambda]qcd/Sqrt[Mpl/(2 t1 t)])^4;
madt[t_] := ma*If[t < t0, m0[t]/ma, 1];


mHDt[t_] := 
  Sqrt[1/2 madt[t]^2 (1 + 
      Rm^2 ma^2/
        madt[t]^2 (1 + Rf^2 + 
         Sqrt[(1 - Rf^2 - 1/Rm^2 madt[t]^2/ma^2)^2 + 4 Rf^2]))];
mLDt[t_] := 
  Sqrt[1/2 madt[t]^2 (1 + 
      Rm^2 ma^2/
        madt[t]^2 (1 + Rf^2 - 
         Sqrt[(1 - Rf^2 - 1/Rm^2 madt[t]^2/ma^2)^2 + 4 Rf^2]))];

    
ti = 1/1000;
tlc = t0 (1 - Rf^2)^(1/4) Sqrt[
  Rm];
tf =  4 tlc; a0 = \[Theta]*fa; phi0 = \[Phi]*fphi;

time = {t, ti, tf};
qvars = {a[t], \[CurlyPhi][t]};

(*DIFFERENTIAL EQUATION*)

eqns = {(a''[t] + 3/(2 t) a'[t]) + t1^2 madt[t]^2 a[t] + 
      t1^2 mphi^2 Rf^2 a[t] + t1^2 mphi^2 Rf phi0/a0 \[CurlyPhi][t] ==
      0, (\[CurlyPhi]''[t] + 3/(2 t) \[CurlyPhi]'[t]) + 
      t1^2 mphi^2 \[CurlyPhi][t] + t1^2 mphi^2 Rf a0/phi0 a[t] == 0, 
    a'[ti] == 0, \[CurlyPhi]'[ti] == 0, 
    a[ti] == 1, \[CurlyPhi][ti] == 1} // Simplify;

solA = NDSolve[eqns, {a, \[CurlyPhi]}, time, StartingStepSize -> 1/20,
     Method -> "BDF", InterpolationOrder -> 2, 
    MaxSteps -> Infinity]; // AbsoluteTiming
```
Improve this question
$\endgroup$
2
  • 2
    $\begingroup$ Actually your code doesn't run. Several parameters are ti be defined! $\endgroup$
    – Ulrich Neumann
    Jan 5 at 10:16
  • $\begingroup$ @UlrichNeumann Thanks a lot for your comment. I tried to copy the above code in my machine and it ran fine, I think it might be the case of the names of variables( curly phi etc) which are not showing up correctly when we copy it to the mathematica notebook. $\endgroup$
    – Mudit Rai
    Jan 6 at 9:57

1 Answer 1

Reset to default
5
$\begingroup$

To simplify code we use substitutions f\[CurlyPhi]->fphi, m\[CurlyPhi]->mphi, ma$t->madt, mH$t->mHDt, mL$t->mLDt, t->t/t1, a->a/a0, \[CurlyPhi]->\[CurlyPhi]/phi0 where a0 = \[Theta] fa; phi0 = \[Phi] fphi. In this variables code can be written as follows

(*VALUES FOR VARIABLES USED*)
fa = 1.13*10^13;
Rf = 0.0104;
\[Theta] = 1;
\[Phi] = 1;
Mpl = 12*10^18;
\[CapitalLambda]qcd = 2/10;
fphi = Rf*fa;


t1 = (Mpl^2/(8 ma \[CapitalLambda]qcd^4))^(1/3);
ma = 57*10^-16 (10^12/fa);

Rm = Rf^4 (3/10 87/10 24/10 1/1000 (fa/10^12)^(3/2))/(11/100);
mphi = Rm ma; t0 = Mpl/(2 \[CapitalLambda]qcd^2)/t1; 
m0[t_] := ma (\[CapitalLambda]qcd/Sqrt[Mpl/(2 t1 t)])^4;

(*Mass matrix eigenvalues*)
madt[t_] := ma If[t < t0, m0[t]/ma, 1];
mHDt[t_] = 
  Sqrt[1/2 madt[t]^2 (1 + 
      Rm^2 ma^2/
        madt[t]^2 (1 + Rf^2 + 
         Sqrt[(1 - Rf^2 - 1/Rm^2 madt[t]^2/ma^2)^2 + 4 Rf^2]))];

mLDt[t_] = 
  Sqrt[1/2 madt[t]^2 (1 + 
      Rm^2 ma^2/
        madt[t]^2 (1 + Rf^2 - 
         Sqrt[(1 - Rf^2 - 1/Rm^2 madt[t]^2/ma^2)^2 + 4 Rf^2]))];

(*Initial and final solving time,increasing the final times makes the \
equation take forever*)
ti = 1/1000;
tf = 250; a0 = \[Theta] fa; phi0 = \[Phi] fphi;
(*DIFFERENTIAL EQUATION*)
eqns = {(a''[t] + 3/(2 t) a'[t]) + t1^2 madt[t]^2 a[t] + 
     t1^2 mphi^2 Rf^2 a[t] + t1^2 mphi^2 Rf phi0/a0 \[CurlyPhi][t] == 
    0, (\[CurlyPhi]''[t] + 3/(2 t) \[CurlyPhi]'[t]) + 
     t1^2 mphi^2 \[CurlyPhi][t] + t1^2 mphi^2 Rf a0/phi0 a[t] == 0, 
   a'[ti] == 0, \[CurlyPhi]'[ti] == 0, 
   a[ti] == 1, \[CurlyPhi][ti] == 1} // Simplify

Now it tak3s about 75s on my laptop

sol = NDSolve[eqns, {a, \[CurlyPhi]}, {t, ti, tf}, 
    MaxSteps -> Infinity]; // AbsoluteTiming

Visualization

Plot[a[t] /. sol[[1]], {t, ti, tf}, PlotRange -> All, Frame -> True, 
 FrameLabel -> {"t", "a"}]

Plot[\[CurlyPhi][t] /. sol[[1]], {t, ti, tf}, PlotRange -> All, 
 FrameLabel -> {"t", "\[CurlyPhi]"}, Frame -> True]

Figure 1

Update 1. Let check equation for a, we have

a[t]*(1.396316111230277*^-13 + 2.0160295703419477*^6*If[t < 53.28920267804457, m0[t]/ma, 1]^2) + 1.396316111230277*^-13*\[CurlyPhi][t] + (3*Derivative[1][a][t])/(2*t) + Derivative[2][a][t] ==0

In this range of parameters we can drop small member 1.396316111230277*^-13*\[CurlyPhi][t], and solve equation for a separately for t<t0 and t>t0 in the form

sol2 = 
 DSolve[{a[t] k t^4 + (3 Derivative[1][a][t])/(
      2 t) + (a^\[Prime]\[Prime])[t] == 0, 
    Derivative[1][a][1/1000] == 0, a[1/1000] == 1}, a, t][[1]]

Out[]= {a -> 
  Function[{t}, -((-Sqrt[k] BesselJ[-(11/12), Sqrt[k]/
          3000000000] BesselJ[-(1/12), (Sqrt[k] t^3)/3] + 
        500000000 BesselJ[-(1/12), (Sqrt[k] t^3)/3] BesselJ[1/12, 
          Sqrt[k]/3000000000] + 
        Sqrt[k] BesselJ[-(13/12), Sqrt[k]/3000000000] BesselJ[1/12, (
          Sqrt[k] t^3)/3] - 
        500000000 BesselJ[-(1/12), Sqrt[k]/3000000000] BesselJ[1/12, (
          Sqrt[k] t^3)/3] - 
        Sqrt[k] BesselJ[1/12, (Sqrt[k] t^3)/3] BesselJ[11/12, Sqrt[k]/
          3000000000] + 
        Sqrt[k] BesselJ[-(1/12), (Sqrt[k] t^3)/3] BesselJ[13/12, Sqrt[
          k]/3000000000])/(10^(3/4) Sqrt[k] t^(
        1/4) (BesselJ[-(11/12), Sqrt[k]/3000000000] BesselJ[-(1/12), 
            Sqrt[k]/3000000000] - 
          BesselJ[-(13/12), Sqrt[k]/3000000000] BesselJ[1/12, Sqrt[k]/
            3000000000] + 
          BesselJ[1/12, Sqrt[k]/3000000000] BesselJ[11/12, Sqrt[k]/
            3000000000] - 
          BesselJ[-(1/12), Sqrt[k]/3000000000] BesselJ[13/12, Sqrt[k]/
            3000000000])))]} 


sol3 = 
 DSolve[a[t] k1 + (3 Derivative[1][a][t])/(
     2 t) + (a^\[Prime]\[Prime])[t] == 0, a, t][[1]]

Out[]= {a -> 
  Function[{t}, (BesselJ[1/4, Sqrt[k1] t] C[1])/t^(1/4) + (
    BesselY[1/4, Sqrt[k1] t] C[2])/t^(1/4)]}

We can compare exact solution sol2 with numerical solution as follows

Show[Plot[a[t] /. sol2 /. k -> 1/4, {t, ti, 10}, PlotRange -> All, 
 PlotPoints -> 200, Frame -> True],Plot[a[t] /. sol[[1]], {t, ti, 10}, PlotStyle -> {Red, Dashed}]]

Figure 2

Update 2. To reduce computation time we can use different strategy at t<t0, and t>t0 as follows

time = {t, ti, t0};
qvars = {a[t], \[CurlyPhi][t]}; 
 solL = NDSolve[eqn, {a, \[CurlyPhi]}, time, StartingStepSize -> 1/20,
    Method -> "Adams", InterpolationOrder -> 2, 
   MaxSteps -> Infinity]; // AbsoluteTiming

It takes about 1s only. Using solL we can compute initial data for the next step

ini = Flatten[{qvars, D[qvars, t]} /. solL /. t -> t0]

eqn1 = {(a''[t] + 3/(2 t) a'[t]) + t1^2 ma^2 a[t] + 
     t1^2 mphi^2 Rf^2 a[t] + t1^2 mphi^2 Rf phi0/a0 \[CurlyPhi][t] == 
    0, (\[CurlyPhi]''[t] + 3/(2 t) \[CurlyPhi]'[t]) + 
     t1^2 mphi^2 \[CurlyPhi][t] + t1^2 mphi^2 Rf a0/phi0 a[t] == 0, 
   a'[t0] == ini[[3]], \[CurlyPhi]'[t0] == ini[[4]], 
   a[t0] == ini[[1]], \[CurlyPhi][t0] == ini[[2]]};

solA = NDSolve[eqn1, {a, \[CurlyPhi]}, {t, t0, tf}, 
    MaxSteps -> Infinity, Method -> "BDF", 
    InterpolationOrder -> 2]; // AbsoluteTiming

It takes about 18 s. We can compare exact solution sol2 with numerical solution solL as follows

Show[Plot[a[t] /. sol2 /. k -> 1/4, {t, ti, 10}, PlotRange -> All, 
 PlotPoints -> 200, Frame -> True],Plot[a[t] /. solL[[1]], {t, ti, 10}, PlotStyle -> {Green, Dashed}]]

Figure 3

Improve this answer
$\endgroup$
12
  • $\begingroup$ Hi Alex, Thanks a lot for your answer. I think it is a good strategy to use dimensionless quantities as implied in your solution, but I think we are still getting O(100) seconds as the time per point. I tried using MaxStepSize, which lead to bogus answer. Using PrecisionGoal->5, it cut down the time by half, still it is a lot as I need to scan over around 10^3 or 10^4 points. Do you have any suggestions which can speed things up considerably? $\endgroup$
    – Mudit Rai
    Jan 6 at 10:00
  • 1
    $\begingroup$ Maybe we can use exact solution to this problem? $\endgroup$
    – Alex Trounev
    Jan 6 at 16:45
  • $\begingroup$ Hi Alex, thanks again for the updates. I think the exact solution obtained by dropping the other term, which resulted in the coupled ODE defeats the purpose of my problem since I would be scanning over possible points where it is very important. $\endgroup$
    – Mudit Rai
    Jan 11 at 10:50
  • $\begingroup$ On another note, the Update 2 part of partitioning the time interval looks promising to me. I wanted to ask if you had defined the t_0 somewhere which you were using in the example, or did you pick some value looking at the function's profile which separates the regions judiciously. Thanks again, these answers are definitely helpful in general for dealing with NDSolve! $\endgroup$
    – Mudit Rai
    Jan 11 at 10:52
  • 1
    $\begingroup$ @MuditRai Your question about Adams and BDF methods is good one, but not suitable for this forum. See, for instance scicomp.stackexchange.com/search?q=Adams+method and scicomp.stackexchange.com/search?q=BDF+method $\endgroup$
    – Alex Trounev
    Feb 3 at 11:43

Your Answer

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

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