Problem with NIntegrate and more general problems with NDSolve

I'm new to Mathematica and I'm receiving little help for a (to me) quite difficult problem I'm trying to solve. I want to solve logarithmic Boltzmann equation for abundance of dark matter. I'm having different problems:

1. I obtain the error "The integrand has evaluated to non-numerical values..." even though I've tried to use NumericQ as read in this forum
2. "Encountered non-numerical value for a derivative at x=0.01" (which I think may be due to the presence of neq as a non evaluated function)

I noted neq isn't evaluated because of the fact I'm defining it with x as a numerical variable, but I need to pass it to NDSolve and I don't know how to do that with this.

Needs["DifferentialEquationsNDSolveUtilities"];
Mx = 100;
gx = 2;
gs = 106.75;
g = 100;
Mp = 10^19;
\[Sigma]v = 10^-10;

Clear[neq, f];
neq[x_?NumericQ, m_] :=
NIntegrate[(gx/(2 Pi Pi)) Exp[-x Sqrt[1 + p p/(m m)]] p p, {p, 0,
Infinity}, MaxRecursion -> 50];
stot[x_] := (2 Pi Pi/45) gs Mx^3 /(x^3);
Yeq[x_, m_] := neq[x, m]/stot[x];
Weq[x_, m_] := Log[Yeq[x, m]];
xmin = 10^-2; xmax = 10^2;
xdom = {x, xmin, xmax};
ODE = {W'[
x] == (0.264 Mp Mx gs \[Sigma]v/(Sqrt[g])) (Exp[
2 Weq[x, m] - W[x]] - Exp[W[x]])/(x^2)};
Y0 = 0;
W0 = -10^10; (*ideally I'd have W0 = -Infinity and Y0 = 0 being W0 = \
log(Y0)*)
BC = {W[xmin] == W0};
Wsol = NDSolve[{ODE, BC}, W[x], xdom, Method -> "StiffnessSwitching"]

Plot[Evaluate[W[x] /. Wsol], {x, 10^-2, 10^2}, PlotRange -> Automatic]


Edit. I'm adding the expected behaviour of the solution from x = 1 to x = xmax

• The first error is due to a wrong definition of neq: the variable m needs to be numerical too, so the definition should be neq[x_?NumericQ, m_?NumericQ] :=.... This however does not solve the all the problems, because another error appears with NDSolve complaining that there are more variable than equations. Commented May 2, 2022 at 14:37
• Also inside the exponential, in Weq there is an undefined m: which constant should it be? Giving a random value to m the program works. In the above comment I forgot to say that in the definition of all of your functions you should add ?NumericQ for all the variables not just in neq. Commented May 2, 2022 at 14:52
• One more thing: Mp and \[Sigma]v differ by 29 orders of magnitude, this is likely to give numerical problems. Commented May 2, 2022 at 14:55
• Hello! Thanks for your observations. I just corrected that little m in Weq(x,m), thanks! I changed my variables to numerical ones as you say, but still it doesn't plot a thing.. does it work to you in that way? Commented May 2, 2022 at 14:57
• Mp and Sigmav are so, the first one is Planck mass and the second one in a specific value of cross section I can't modify. The stiffness of the equation is part of the problem I'm trying to overcome Commented May 2, 2022 at 14:58

Here is the improved code. I have changed the method of NDSolveValue to Method -> {"StiffnessSwitching", Method -> {"ExplicitRungeKutta", Automatic}} and increased the AccuracyGoal. Better results require higher accuracy, but will be much more time consuming (The code below required about a minute on my computer). So here is the code and the plot.

Needs["DifferentialEquationsNDSolveUtilities"];
Mx = 100;
gx = 2;
gs = 106.75;
g = 100;
Mp = 10^19;
\[Sigma]v = 10^-10;

Clear[neq];
neq[x_?NumericQ, m_?NumericQ] :=
NIntegrate[(gx/(2 Pi Pi)) Exp[-x Sqrt[1 + p p/(m m)]] p p, {p, 0,
Infinity}];
stot[x_?NumericQ] := (2 Pi Pi/45) gs Mx^3/(x^3);
Yeq[x_?NumericQ, m_?NumericQ] := neq[x, m]/stot[x];
xmin = 10^-1; xmax = 10^2;
xdom = {x, xmin, xmax};
ODE = {Y'[
x] == (0.264 Mp Mx gs \[Sigma]v/(Sqrt[g])) (Yeq[x, Mx]^2 -
Y[x]^2)/(x^2)};
BC = {Y[xmin] == 10^-6};
Wsol = NDSolveValue[{ODE, BC}, Y, xdom,
Method -> {"StiffnessSwitching",
Method -> {"ExplicitRungeKutta", Automatic}}, AccuracyGoal -> 20]
LogLogPlot[Wsol[x], {x, 0.1, 100}, PlotRange -> All]


Plot:

Ok here is a possible solution. Since W0 is too high, you could write the ODE in terms of Y[x]=Exp[W[x]] setting a small number for the initial condition. For example you could do something like this:

Needs["DifferentialEquationsNDSolveUtilities"];
Mx = 100;
gx = 2;
gs = 106.75;
g = 100;
Mp = 10^19;
\[Sigma]v = 10^-10;

Clear[neq, f];
neq[x_?NumericQ, m_?NumericQ] :=
NIntegrate[(gx/(2 Pi Pi)) Exp[-x Sqrt[1 + p p/(m m)]] p p, {p, 0,
Infinity}, MaxRecursion -> 50];
stot[x_?NumericQ] := (2 Pi Pi/45) gs Mx^3/(x^3);
Yeq[x_?NumericQ, m_?NumericQ] := neq[x, m]/stot[x];
xmin = 10^-2; xmax = 10^2;
xdom = {x, xmin, xmax};
ODE = {Y'[x] ==
Y[x] (0.264 Mp Mx gs \[Sigma]v/(Sqrt[g])) (Yeq[x, Mx]/Y[x] -
Y[x])/(x^2)};
BC = {Y[xmin] == 0.1};
Wsol = NDSolveValue[{ODE, BC}, Y[x], xdom, Method -> "StiffnessSwitching"]

Plot[Wsol,xdom]


The plot is ugly:

But notice that in the NDSolveValue code I set m=Mx just to put in a value, since you did not provide the correct constant, so the result may vary. Also I set Y[xmin] to 0.1, probably you should play a bit with this boundary.

• Hello! First of all, thank you so much for helping me with this problem, it means a lot. Now m = Mx is correct as you guessed, but the boundary condition needs to be given for Y(xmin) = 0 or at least very close to it. I was working with W(x) = log(Y(x)) to make the integration more efficient since Y spans over more than 40 orders of magnitudes. Using the BC I told you I obtain something nosense given that the expected behaviour of Y(x) should be decreasing. I understand this is a particularly hard problem to deal with. I'm adding the (expected) plot we should obtain from x = 1 to x = xmax Commented May 2, 2022 at 16:11
• @Fredrigo6 See my update. I will try later to see if can do better Commented May 2, 2022 at 17:05
• Hello! Thanks a lot for your help. I'm going to bed now so I'll talk to you tomorrow! Commented May 2, 2022 at 17:47
• @Fredrigo6. See the update. Probably rescaling the parameters would require smaller accuracy, but if you can't change them this is a solution. Commented May 2, 2022 at 18:03
• This kinda works! I guess it's just a difficult problem. Imagine trying to solve it from x = 0.001 and with boundary condition Y[xi] = 0. thanks for you help though! I'll accept the answer, but if you are intrigued by it I'll continue working on it so in case let's continue the conversation! (How can I accept your answer?) Commented May 3, 2022 at 6:49