3
$\begingroup$

I have the code below. If you look at the plot, you can see a strange behaviour. The function goes to zero as λ->0 but you can see the plateau below approximately 5e-8. Yet, if I evaluate the function at 1e-8 I get basically zero. What is the source of this behaviour?

enter image description here

ch = QuantityMagnitude[UnitConvert[Quantity[1, "PlanckConstant"], "SIBase"]];
cc = QuantityMagnitude[UnitConvert[Quantity[1, "SpeedOfLight"], "SIBase"]];
ckB = QuantityMagnitude[UnitConvert[Quantity[1, "BoltzmannConstant"], "SIBase"]];
TSun = 5800;(* Kelvins *)
TEarth = 255;(* Kelvins *)

Bλ[T_, λ_] := (2 ch*cc^2)/λ^5 1/(Exp[ch/(ckB*T) cc/λ] - 1);

λH = 1000*10^-6;(* "Meters" *)
λL = 10.0*10^-9;(* "Meters" *)

LogLogPlot[Bλ[TEarth, λ], {λ, λL, λH}, Frame -> True]

Bλ[TEarth, 1.0*10^-8]
Improve this question
$\endgroup$
1
  • 2
    $\begingroup$ The workaround of choosing the properly scaled units here (e.g. nanometers) should alleviate this problem. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Feb 15, 2016 at 2:11

2 Answers 2

Reset to default
3
$\begingroup$

I think it is because $10^{-2427}$ is beyond the $MinMachineNumber, $2.2250738585072014 \times 10^{-308}$.

LogLogPlot[
 Bλ[TEarth, λ], {λ, λL, λH}, 
 Frame -> True, PlotRange -> All]

enter image description here

Improve this answer
$\endgroup$
3
$\begingroup$

I don't know if you can call this a bug, but it seems buggy to me. Let's take a look at the data that is actually being plotted, which you can extract using a Reap and Sow combination (thanks to @user21 for showing me this trick)

Reap[LogLogPlot[Sow[{λ, Bλ[TEarth, λ]}]; 
     Bλ[TEarth, λ], {λ, λL, λH}, 
     Frame -> True]][[2, 1, ;; 20]] // TableForm

enter image description here

That would seem to be the big glaring problem, right there, that point of {x,f[x]} = {-18.4207, 1.83338*10^-17}. But I'm not actually convinced, I think that is a red herring, because you get the same x value if you replace Bλ[TEarth, λ] with λ^2.

About the only thing we do learn from that table is that Mathematica doesn't seem to have a problem with really small numbers.

This is apparently one of those cases where you need to use Evaluate in the argument to Plot,

GraphicsRow[{
  LogLogPlot[
   Bλ[TEarth, λ], {λ, λL, λH},
    Frame -> True],
  LogLogPlot[
   Evaluate@
    Bλ[
     TEarth, λ], {λ, λL, λH}, 
   Frame -> True]
  }, ImageSize -> 700]

enter image description here

This is a simple workaround that you can use all the time without drawback.

Working Precision giving crazy results

It was suggested to check out the results of using WorkingPrecision on the plot, so look at this mess:

LogLogPlot[
   Bλ[TEarth, λ], {λ, λL, λH},
    Frame -> True, PlotRange -> All, WorkingPrecision -> #, 
   ImageSize -> 300] & /@ {MachinePrecision/2, MachinePrecision, 
  2 MachinePrecision, $MachinePrecision}

enter image description here

The only setting that was even close to correct was the default setting. It seems the only real answer here is to wrap the function with Evaluate before plotting.

Edit Here is my original answer below, which I'm leaving in only because it shows off the version 10 function PlanckRadiationLaw, which is a little slow for sure.

@J.M. made the point of using nanometers instead of meters, and my first thought was to use atomic units instead of SI units, but you still run into the same problem when taking a LogLog plot.

Your code is just fine, but I thought I'd show off a version 10 function here,

LogLogPlot[{
  QuantityMagnitude[
   PlanckRadiationLaw[Quantity[5778, "Kelvins"], 
    Quantity[x, "Nanometers"]]],
  QuantityMagnitude[
   PlanckRadiationLaw[Quantity[255, "Kelvins"], 
    Quantity[x, "Nanometers"]]]}
 , {x, 10, 10^6}, 
 AxesLabel -> {Quantity[None, "Nanometers"], 
   Quantity[None, "Watts"/("Hertz"*"Meters"^2*"Steradians")]}]

enter image description here

shows the same problem, even with proper units. But if you adjust the lowest wavelength to 100 nm, then you get this:

enter image description here

Improve this answer
$\endgroup$
4
  • $\begingroup$ I understand how to correct it, I am more curious why the behaviour is like this, ie. why Mathematica shows this apparently wrong limit for the function? Also, notice that Mathematica does not issue any warning that the number is too small and it is confusing, I think there can be more complicated functions where the result is not so apparently wrong and you can produce nonsense without noticing it. Also, I am curious, why this really happens, is the issue with comparing the numbers that are beyond the $MinMachineNumber or is it some type of overflow? $\endgroup$
    – atapaka
    Commented Feb 15, 2016 at 14:08
  • $\begingroup$ @leosenko - you may have to write to the developers to get a good answer. You can rephrase the question to leave out the extraneous stuff, like the units and planck's constant. You can just define a function as f[x_] := 10^-16 x^-5/(Exp[5 10^-5/x]); and then compare the results of LogLogPlot[f[x], {x, 10^-8, 10^-3}] to LogLogPlot[10^-16 x^-5/(Exp[5 10^-5/x]), {x, 10^-8, 10^-3}] That would make it a minimal example that demonstrates the problem. $\endgroup$
    – Jason B.
    Commented Feb 15, 2016 at 15:13
  • $\begingroup$ Changing the setting of WorkingPrecision is also useful: LogLogPlot[Bλ[TEarth, λ], {λ, λL, λH}, Frame -> True, PlotRange -> All, WorkingPrecision -> $MachinePrecision]. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Feb 15, 2016 at 17:20
  • $\begingroup$ @J.M. - That's crazy, just makes me even more confused. See the edit $\endgroup$
    – Jason B.
    Commented Feb 16, 2016 at 13:06

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.