4
$\begingroup$

This is my first post on this site, please forgive me if I screw up to copy my code in a nice way..

I am trying to solve a set of 3 coupled differential equations with NDSolve, where one of the terms in the equations is obtained by NIntegrate. The set of equations looks as follows

NDSolve[{D[X[A], A] == (A^(7/2)/TRHb^5)*((Sqrt[3]*Mpl*RF[A,10^8,10^12])/Sqrt[Φ[A] + R[A]/A + (X[A]*Sqrt[mDM^2 + 9*κb^2*(A^(-3/2))^(1/2)])/TRHb]),

         D[Φ[A], A] == (-(A^(1/2)/TRHb^2))*((Sqrt[3]*Mpl*ΓbSM)/Sqrt[Φ[A] + R[A]/A + (X[A]*Sqrt[mDM^2 + 9*κb^2*(A^(-3/2))^(1/2)])/TRHb]), 

         D[R[A], A] == (ΓbSM*Φ[A] - (1/TRHb^4)*(2*Sqrt[mDM^2 + 9*κb^2*(A^(-3/2))^(1/2)]*RF[A,10^8,10^12]*A^3))*((Sqrt[3]*Mpl*A^(3/2))/(TRHb^2*Sqrt[Φ[A] + R[A]/A + (X[A]*Sqrt[mDM^2 + 9*κb^2*(A^(-3/2))^(1/2)])/TRHb])),

        X[10^(-10)] == 0,
        Φ[10^(-10)] == ΦI,
        R[10^(-10)] == 0},

        {X, Φ, R}, {A, 10^(-10), 10^4}]}

The function RF is the one which I need to solve with NIntegrate. The final goal is to solve this set of equations for X[A,mDM,λ]. In order to do so, I am trying to solve these equations on a "grid" of different mDM and λ values. Afterwards I fix A and use Interpolation in order to get a nice smooth solution in mDM and λ. But I already struggle to get a solution just for a single combination of mDM and λ (mDM = 10^8, λ = 10^12). I guess the main problem is in the function RF. It is defined as follows

RF[(A_)?NumericQ, (mDM_)?NumericQ, (mt_)?NumericQ, (αG_)?NumericQ, (λ_)?NumericQ]:=

 (4*Pi/(2048*Pi^6))*NIntegrate[

Norm[(Sqrt[1-(4*mDM^2)/s]*((mDM^2*s^3*(1-(4*mDM^2)/s))*(4 + 12*αG^2)^2)*κb*(A^(-3/2))^(1/4)*Sqrt[s]*BesselK[1, Sqrt[s]/(κb*(A^(-3/2))^(1/4))])/(λ^4*((s - mt^2)^2 + mt^2*((mt^3*((1/8)*Sqrt[1 - (4*89^2)/mt^2] + (3*αG^2)/16))/(Pi*λ^2) + (mt^3*Sqrt[1 - (4*mDM^2)/mt^2]*mDM^2*(1 - (4*mDM^2)/mt^2))/((Pi*λ^2)*8*mt^2))^2))],

                    {s, 4*mDM^2, Infinity}, Exclusions -> {mt^2, ma^2}]

As a function of A, RF[A] should look like FIG1. in [ https://arxiv.org/abs/1811.01947 ]. But I fail to reproduce this form in the rage mDM = [10^8 - 10^14]. As soon as NIntegrate hits the point s ~ mt^2, I get problems. I was trying to play around with the NIntegrate settings, but nothing really helped. I obtained the best solutions with "DoubleExponential" and "Trapezoidal", but in both cases I get Singaluarity errors. I guess I need to get first an error free solution in NIntegrate before I can use RF in NDSolve. How can I get a nice integration for the entire range of mDM?

The parameters I used were

ma = 0;αS = 1; αV = 1/Log[10^5]; αH = 1; αG = 1/Log[10^5]; βG = 0; βV = 0; αχ = 1; αχV = 1; βχ = 0; βχ5 = 0; g = 100; V = 10^5; Mpl = 2*10^18; κinf = 10^17; κb = 9*10^8; ΦI = 4*10^43; ΓinfSM = 4*10^8; ΓbSM = 7/10^11; minf = 6*10^14; mt = 2*10^10; TRHinf = 2*10^13; TRHb = 7*10^4;

Update

Here are some plots of the function RF[A] for the parameters mDM = 10^8, 10^9, 10^10, 10^11, 10^12, 10^13. Inside NIntegrate I used

{s, 4 mDM^2, mt^2, Infinity}, Method -> "DoubleExponential", WorkingPrecision -> 10, AccuracyGoal -> 10.

Timing[Plot[RF[A, mDM, mt, αG, λ], {A, 10^(-10), 10^10}, ScalingFunctions -> {"Log", "Log"}, PlotRange -> Automatic, PerformanceGoal -> "Quality"]]

enter image description here

enter image description here

enter image description here

enter image description here

enter image description here

enter image description here

And I get the following error for the case of mDM=10^12

NIntegrate::ncvi: NIntegrate failed to converge to prescribed accuracy after 9 iterated refinements in s in the region {{4.00000000000000000000000000000000000000000000000000000000000*10^20,\[Infinity]}}. NIntegrate obtained 9.2531746437536964660643301435717812896328729569142552371394`60.*^54 and 2.60922345176519777170004482614311324727140650854448718254727`60.*^53 for the integral and error estimates.

I guess before I use RF[A] as an input for NDSolve I need to get it as a smooth function for the entire parameter space.

Improve this question
$\endgroup$
20
  • $\begingroup$ I recommend stripping out the unimportant parts and simplifying the problem until you have the kernel of it. It's awkward picking through your code at the moment. $\endgroup$
    – MikeY
    Mar 4, 2020 at 17:56
  • $\begingroup$ I just cleared up the code a little bit, I hope it is easier to understand now. $\endgroup$
    – ikarus
    Mar 4, 2020 at 19:07
  • $\begingroup$ @ikarus What do you want to calculate here? $\endgroup$
    – Alex Trounev
    Mar 4, 2020 at 20:27
  • $\begingroup$ @AlexTrounev I am trying to solve the 3 coupled Boltzmann equations, which describe the production of Dark Matter during Inflation. X is dark matter, Phi the inflaton and R is radiation. The function RF is the production rate of DM, i.e basically the decay of radiation to dark matter $\endgroup$
    – ikarus
    Mar 4, 2020 at 20:45
  • $\begingroup$ @ikarus This is clear from the article. What is your problem? I checked the code on version 12 for Windows 10. Everything works. $\endgroup$
    – Alex Trounev
    Mar 4, 2020 at 20:50

1 Answer 1

Reset to default
2
$\begingroup$

I can offer one working code. I checked this for a parameter set, but not sure if this works for any parameter sets

ma = 0; \[Alpha]S = 1; \[Alpha]V = 
 1/Log[10^5]; \[Alpha]H = 1; \[Alpha]G = 
 1/Log[10^5]; \[Beta]G = 0; \[Beta]V = 0; \[Alpha]\[Chi] = 1; \
\[Alpha]\[Chi]V = 1; \[Beta]\[Chi] = 0; \[Beta]\[Chi]5 = 0; g = 100; \
V = 10^5; Mpl = 2*10^18; \[Kappa]inf = 10^17; \[Kappa]b = 
 9*10^8; \[CapitalPhi]I = 4*10^43; \[CapitalGamma]infSM = 
 4*10^8; \[CapitalGamma]bSM = 7/10^11; minf = 6*10^14; mt = 
 2*10^10; TRHinf = 2*10^13; TRHb = 7*10^4; 
RF[(A_)?NumericQ, (mDM_)?NumericQ, (mt_)?NumericQ, (ma_)?
   NumericQ, (\[Alpha]H_)?NumericQ, (\[Alpha]G_)?
   NumericQ, (\[Alpha]\[Chi]_)?NumericQ, (\[Alpha]\[Chi]V_)?
   NumericQ, (\[Beta]\[Chi]_)?NumericQ, (\[Beta]\[Chi]5_)?
   NumericQ, (\[Beta]G_)?NumericQ, (\[Lambda]_)?
   NumericQ] := ((4*Pi)*
    NIntegrate[
     Norm[(Sqrt[
          1 - (4*mDM^2)/
            s]*(\[Alpha]\[Chi]^2*(mDM^2*
             s^3*(1 - (4*mDM^2)/s))*(4*\[Alpha]H^2 + 
              12*\[Alpha]G^2)^2)*\[Kappa]b*(A^(-3/2))^(1/4)*Sqrt[s]*
         BesselK[1, 
          Sqrt[s]/(\[Kappa]b*(A^(-3/2))^(1/4))])/(\[Lambda]^4*((s - 
              mt^2)^2 + 
           mt^2*((mt^3*((1/8)*\[Alpha]H^2*
                    Sqrt[1 - (4*89^2)/mt^2] + (3*\[Alpha]G^2)/
                    16))/(Pi*\[Lambda]^2) + (mt^3*
                  Sqrt[1 - (4*mDM^2)/mt^2]*\[Alpha]\[Chi]V^2*
                  mDM^2*(1 - (4*mDM^2)/mt^2))/((Pi*\[Lambda]^2)*8*
                  mt^2))^2))], {s, 4*mDM^2, Infinity}, 
     Exclusions -> {mt^2, ma^2}, AccuracyGoal -> 10])/(2048*Pi^6);
x[m_, l_] := 
 Module[{mDM = m, lambda = l}, 
  rf[A_] := 
   RF[A, mDM, mt, 
    ma, \[Alpha]H, \[Alpha]G, \[Alpha]\[Chi], \[Alpha]\[Chi]V, \
\[Beta]\[Chi], \[Beta]\[Chi]5, \[Beta]G, lambda];
  eq = {D[X[A], 
      A] == (A^(7/2)/TRHb^5)*((Sqrt[3]*Mpl*rf[A])/
        Sqrt[\[CapitalPhi][A] + 
          R[A]/A + (X[A]*
             Sqrt[mDM^2 + 9*\[Kappa]b^2*(A^(-3/2))^(1/2)])/TRHb]), 
    D[\[CapitalPhi][A], 
      A] == (-(A^(1/2)/TRHb^2))*((Sqrt[3]*Mpl*\[CapitalGamma]bSM)/
        Sqrt[\[CapitalPhi][A] + 
          R[A]/A + (X[A]*
             Sqrt[mDM^2 + 9*\[Kappa]b^2*(A^(-3/2))^(1/2)])/TRHb]), 
    D[R[A], A] == (\[CapitalGamma]bSM*\[CapitalPhi][
          A] - (1/TRHb^4)*(2*
           Sqrt[mDM^2 + 9*\[Kappa]b^2*(A^(-3/2))^(1/2)]*rf[A]*
           A^3))*((Sqrt[3]*Mpl*A^(3/2))/(TRHb^2*
          Sqrt[\[CapitalPhi][A] + 
            R[A]/A + (X[A]*
               Sqrt[mDM^2 + 9*\[Kappa]b^2*(A^(-3/2))^(1/2)])/TRHb]))};
  ic = {X[10^(-10)] == 0, \[CapitalPhi][10^(-10)] == \[CapitalPhi]I, 
    R[10^(-10)] == 0};
  f = NDSolveValue[{eq, ic}, X[A], {A, 10^(-10), 10^4}]; f]
    {Plot[Evaluate[x[10^8, 10^13]], {A, 10^(-10), 10^4}, 
  AxesLabel -> {"A", "X"}], 
 Plot[Evaluate[x[10^8, 10^12]], {A, 10^(-10), 10^4}, 
  AxesLabel -> {"A", "X"}], 
 Plot[Evaluate[x[10^9, 10^12]], {A, 10^(-10), 10^4}, 
  AxesLabel -> {"A", "X"}]}

Figure 1

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thanks a lot for the code, it works indeed! But in order to reproduce Fig1. in the reference, I need to use Method -> "DoubleExponential" inside NIntegrate. Otherwise I do not get the right slope in the low A (or equivalently high T) regime of RF[A]. But when I use "DoubleExponential" I still get convergence/singularity errors. How can I get rid of those errros? $\endgroup$
    – ikarus
    Mar 9, 2020 at 11:15
  • $\begingroup$ To reproduce Fig. 1, we must obtain an exact expression of the integral, or at least a piecewise function.Plot and NDSolve use RF different ways. $\endgroup$
    – Alex Trounev
    Mar 9, 2020 at 11:28
  • $\begingroup$ Aah, I see. So you mean even though Plot[RF[A]] looks different then in Fig1. I can still trust it inside NDSolve? $\endgroup$
    – ikarus
    Mar 9, 2020 at 11:46

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.