0
$\begingroup$

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["DifferentialEquations`NDSolveUtilities`"];
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]

Thanks for your time!

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

$\endgroup$
7
  • $\begingroup$ 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. $\endgroup$
    – mattiav27
    May 2, 2022 at 14:37
  • $\begingroup$ 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. $\endgroup$
    – mattiav27
    May 2, 2022 at 14:52
  • $\begingroup$ One more thing: Mp and \[Sigma]v differ by 29 orders of magnitude, this is likely to give numerical problems. $\endgroup$
    – mattiav27
    May 2, 2022 at 14:55
  • $\begingroup$ 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? $\endgroup$
    – Fredrigo6
    May 2, 2022 at 14:57
  • $\begingroup$ 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 $\endgroup$
    – Fredrigo6
    May 2, 2022 at 14:58

1 Answer 1

3
$\begingroup$

FINAL ANSWER

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["DifferentialEquations`NDSolveUtilities`"];
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:

enter image description here

OLD ANSWER

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["DifferentialEquations`NDSolveUtilities`"];
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:

enter image description here

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.

$\endgroup$
8
  • $\begingroup$ 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 $\endgroup$
    – Fredrigo6
    May 2, 2022 at 16:11
  • $\begingroup$ @Fredrigo6 See my update. I will try later to see if can do better $\endgroup$
    – mattiav27
    May 2, 2022 at 17:05
  • $\begingroup$ Hello! Thanks a lot for your help. I'm going to bed now so I'll talk to you tomorrow! $\endgroup$
    – Fredrigo6
    May 2, 2022 at 17:47
  • $\begingroup$ @Fredrigo6. See the update. Probably rescaling the parameters would require smaller accuracy, but if you can't change them this is a solution. $\endgroup$
    – mattiav27
    May 2, 2022 at 18:03
  • $\begingroup$ 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?) $\endgroup$
    – Fredrigo6
    May 3, 2022 at 6:49

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.