2

I'm trying to plot two curves but I'm having issues with the lower end of the range:

Clear[t, a0, \[Alpha], \[Delta], \[CapitalTheta], lmin, lmin2, ebar, \
nmin, nmin2]
a0 = 2.4*10^-10;
\[Alpha] = 10;
\[Delta] = 1/100;
\[CapitalTheta] = 0.1;
lmin = a0*nmin;
lmin2 = a0*nmin2;
ebar = ((t/\[CapitalTheta])^2)*
   Integrate[x/(Exp[x] - 1), {x, 0, \[CapitalTheta]/t}];
nmin = (4*\[CapitalTheta]*\[Alpha])/(t*
      ebar)*((ebar/\[Alpha]) + (1/4))^2;
nmin2 = (2 \[Alpha]/\[Delta])*(\[CapitalTheta]/t)*ebar;
LogLogPlot[{nmin, nmin2}, {t, 10^-4, 10^3}, PlotRange -> All]

I keep getting one error and the plot won't show the values for t<10^-2. How can I fix that?

Rodrigo
  • 1,482
  • 9
  • 13

2 Answers2

4

Instead of using Integrate use NIntegrate, since you want numerical solutions. Also, there is a variable of t in integration, for that, you can use := (Set Delayed).

Clear[t, a0, \[Alpha], \[Delta], \[CapitalTheta], lmin, lmin2, ebar, nmin, nmin2]
a0 = 2.4*10^-10;
\[Alpha] = 10;
\[Delta] = 1/100;
\[CapitalTheta] = 0.1;
lmin = a0*nmin;
lmin2 = a0*nmin2;
ebar[t_] := ((t/\[CapitalTheta])^2)*
NIntegrate[x/(Exp[x] - 1), {x, 0, \[CapitalTheta]/t}];
nmin[t_] := (4*\[CapitalTheta]*\[Alpha])/(t*ebar[t])*((ebar[t]/\[Alpha]) + (1/4))^2;
nmin2[t_] := (2 \[Alpha]/\[Delta])*(\[CapitalTheta]/t)*ebar[t];
LogLogPlot[{nmin[t], nmin2[t]}, {t, 10^-4, 10^3}, PlotRange -> All]

enter image description here

sslucifer
  • 541
  • 2
  • 13
4

The underflow is not a significant error. There is a numerics issue with $\log(1+x)$ when $|x| <\mkern-4mu< 1$. Since $\log(1+x)$ is asymptotic to $x$ as $x\rightarrow0$, underflow in $x$ corresponds to underflow in $\log(1+x)$ so there's little to be done. If it's a problem, then use a higher WorkingPrecision, which just pushes back the underflow problem (hopefully to a point at which it is negligible).

The real problem that shows up in the plot is loss of precision in the machine-precision argument 1 + x. Since floating-point has a fixed-width mantissa, when x is smaller than 1, significant bits are lost; when x is much smaller than 1, many signficant bits are lost. Since the result of Log[1 + x] is approximately x, that loss of significance can be, well, significant. There is function log1p in many math libraries that handles this. It's available in Mathematica as Internal`Log1p[], but it is not in the System` context; see Elegant high precision `log1p`?

We can turn off the warning message for a session with Off[] (until turned back on with On[]) or temporarily for a computation with Quiet[].

Clear[t, a0, α, δ, Θ, lmin, lmin2, ebar, nmin, nmin2]
a0 = 2.4*10^-10;
α = 10;
δ = 1/100;
ebar = ((t/Θ)^2)*
    Integrate[x/(Exp[x] - 1), {x, 0, Θ/t}, 
     Assumptions -> Θ > 0 && t > 0] /. 
   Log[x_] :> Internal`Log1p[x - 1 // Simplify];
Θ = 0.1;
nmin = (4*Θ*α)/(t*ebar)*((ebar/α) + (1/4))^2;
nmin2 = (2 α/δ)*(Θ/t)*ebar;
Quiet[
 LogLogPlot[{nmin, nmin2}, {t, 10^-4, 10^3}, PlotRange -> All],
 General::munfl]

enter image description here

Goofy
  • 2,767
  • 10
  • 12