1
$\begingroup$

I have the following program to solve for the neutron to baryon ratio during BBN. I've replaced many of the derived values with constants in order to focus on the problem:

ClearAll[rate, hubble, equation, initialConditions, fraction]
rate[x_] := (0.28*(12 + 6*x + x^2))/x^5
hubble[x_] := 1.73*^-13/Sqrt[x]
equation := Derivative[1][fraction][x] == (rate[x]/(x*hubble[x]))*((1 - fraction[x])/E^x - fraction[x]); 
initialConditions := fraction[0.12] == 0.46
fraction[x_] = fraction[x] /. NDSolve[{equation, initialConditions}, fraction[x], {x, 0.3, 100}]
LogLogPlot[2*fraction[x], {x, 0.65, 65}]

In the domain where $x\approx1$, this ODE goes berserk with these error messages:

NDSolve::ndsz: At x == 0.12, step size is effectively zero; singularity or stiff system suspected.

When I try to plot the ODE solution, I get:

$IterationLimit::itlim: Iteration limit of 4096 exceeded.

Here is an image of the progress of the ODE solution as a function of hubble. I would expect the solution to the ODE to have the shape of the fifth plot in this series as hubble approaches zero.

enter image description here

I've no idea what this means or how to fix it. It appears to be an issue with precision. How can I fix this? My target is to get this section of code working when $hubble=10^{-13}$.

$\endgroup$
6
  • $\begingroup$ Can you run to just before the problem and then plot the solution. Sounds like you have hit a singularity. Look to see if the solution is going off to infinity. $\endgroup$ – Hugh Mar 28 at 15:02
  • $\begingroup$ I'm afraid I don't understand what you are asking. The problem occurs during the evaluation of the ODE. I've got nothing to plot before I have the solution to the ODE. $\endgroup$ – Quarkly Mar 28 at 15:06
  • $\begingroup$ Try reducing the upper limit of the solution range to about 0.12 and plot the solution to find the behaviour as it approaches 0.12 $\endgroup$ – mikado Mar 28 at 15:37
  • $\begingroup$ The function fraction approaches 0.5 as x goes to zero. $\endgroup$ – Quarkly Mar 28 at 15:49
  • $\begingroup$ MMA 12.1 If I run your code, I get a nice solution, no error. Question: is hubble[x_] := 10^(-4)o.k.? $\endgroup$ – Daniel Huber Mar 28 at 15:58
2
$\begingroup$

The computation in the question can be performed by using a higher WorkingPrecison, which in turn requires rationalizing all decimals and employing the option, Method -> "StiffnessSwitching", to avoid the inexplicable NDSolveValue::nderr (Error test failure).

rate[x_] := Rationalize[(0.28*(12 + 6*x + x^2))/x^5, 0]
hubble[x_] := Rationalize[1.73*^-13/Sqrt[x], 0]
equation := Derivative[1][fraction][x] == (rate[x]/(x*hubble[x]))*
    ((1 - fraction[x])/E^x - fraction[x]);
initialConditions := fraction[12/100] == 46/100
sn = NDSolveValue[{equation, initialConditions}, fraction[x], {x, 12/100, 100}, 
    WorkingPrecision -> 45, Method -> "StiffnessSwitching"]
LogLogPlot[2*sn, {x, 0.65, 65}, PlotRange -> {.8 10^-3, 1}, 
    ImageSize -> Large, LabelStyle -> {15, Bold, Black}]

enter image description here

Note that WorkingPrecision -> 30 may appear to give a reasonable plot, but produces inaccurate results when integrating over {x, 12/100, 3/10}, which NDSolve must do to apply the boundary condition. Note also that fraction[x] has an Accuracy of only about -15, and improving it by using WorkingPrecision -> 60 is interminably slow.

$\endgroup$
3
  • $\begingroup$ what purpose does using fractions like $\frac{12}{100}$ instead of real numbers serve? Those constants I placed in the original post (in order to simplify the issue) are actually the results of other, non-trivial functions. Do I have to 'rationalize' them as well? $\endgroup$ – Quarkly Mar 29 at 19:20
  • $\begingroup$ @Quarkly, NDSolve with WorkingPrecision set to a larger value than $MachinePrecision refuses to run, if the the precisions of numbers given to NDSolve are not at least as large as the value of WorkingPrecision. I use Rationalize, because it is easy and gives infinite precision. By the way, I may have a different approach, but I need to think about it a bit. $\endgroup$ – bbgodfrey Mar 29 at 19:32
  • $\begingroup$ This looks very promising. Thank you for the insight into 'Rationalize'. I had tried it, but didn't understand fully how it worked, so I didn't follow it through. $\endgroup$ – Quarkly Mar 29 at 21:40
2
$\begingroup$

Here a version whithout multiple redefinitions of fraction[]:

rate[x_] := ( 28/100*(12 + 6*x + x^2))/x^5
hubble[x_] := hub
equation :=Derivative[1][fraction][
x] == (rate[x]/(x*hubble[x]))*((1 - fraction[x])/E^x - fraction[x]);

initialConditions := fraction[0.12] == 46/100

 fractionN =ParametricNDSolveValue[{equation, initialConditions}, fraction , {x, 0.3, 100}, hub, Method -> "StiffnessSwitching",WorkingPrecision -> 25]

 LogLogPlot[Table[2*fractionN[hub][x], {hub, Table[10^-n, {n, 1, 6}]}] //Evaluate, {x, 0.3, 100}, PlotLegends -> N[Table[10^-n, {n, 1, 6}]]]//Quiet

enter image description here

$\endgroup$
8
  • $\begingroup$ Does the solution for $hubble=10^{-6}$ look at all physical? I appreciate the effort, but that looks to me like a function on the verge of a nervous breakdown. $\endgroup$ – Quarkly Mar 28 at 17:20
  • $\begingroup$ Yes , see my modified answer(WorkingPrecision->25 ) $\endgroup$ – Ulrich Neumann Mar 28 at 17:26
  • $\begingroup$ Sorry, but after researching the ParametricNDSolveValue function, this doesn't do what I need it to do. While the function hubble returns a value on the order of $10^{-13}$, the function, as well as the rate function are both dependent on x and can't be solved outside of the ODE. Trying to force this answer back into NDSolve just brings me back to the original problem. $\endgroup$ – Quarkly Mar 28 at 20:57
  • $\begingroup$ @Quarkly Your proposed edit in the review queue is hard to read (not your fault, it's just formatted in a narrow column). It looked like you were suggesting a change to one line of code in effect. You could do that in a comment, instead of adding such a long edit to this post. (Or I misread it, which is the reason for my first remark.) $\endgroup$ – Michael E2 Mar 28 at 22:20
  • 1
    $\begingroup$ Well, Ulrich can decide. He can review the edit. $\endgroup$ – Michael E2 Mar 29 at 0:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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