# Problem with NIntegrate inside NDSolve

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


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.253174643753696466064330143571781289632872956914255237139460.*^54 and 2.6092234517651977717000448261431132472714065085444871825472760.*^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.

• 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. – MikeY Mar 4 '20 at 17:56
• I just cleared up the code a little bit, I hope it is easier to understand now. – ikarus Mar 4 '20 at 19:07
• @ikarus What do you want to calculate here? – Alex Trounev Mar 4 '20 at 20:27
• @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 – ikarus Mar 4 '20 at 20:45
• @ikarus This is clear from the article. What is your problem? I checked the code on version 12 for Windows 10. Everything works. – Alex Trounev Mar 4 '20 at 20:50

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


• 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? – ikarus Mar 9 '20 at 11:15
• 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. – Alex Trounev Mar 9 '20 at 11:28
• Aah, I see. So you mean even though Plot[RF[A]] looks different then in Fig1. I can still trust it inside NDSolve? – ikarus Mar 9 '20 at 11:46