Infinite expression 1/0. encountered - caused by precision?

Question

Planck's law dependent on frequency rho is as follows:

   B[T_, h_, rho_, k_, c_] := (2 h rho^3)/c^2 1/(Exp[h rho/(T k)] - 1)

As you can see, the denominator can only be zero for rho = 0. All the rest are nonzero constants. Curiously if we try to plot it

Plot[B[5800, 6.626 10^-34, rho, 1.38 10^-23, 299792458], {rho, 
  0.1 10^-6, 3 10^-6}]

It will return "Infinite expression 1/0. encountered" and the plot will be blank. In fact, substituting any value for rho appears to result in this error (out = ComplexInfinity). I have no experience with precision, but I don't see any other possible reason. Anyway to solve this?

Answer 1

To understand what is going on, look at this:

N[b[5800, 6626 10^-37, rho, 138 10^-25, 299792458]]

look at the exponent in the second term in the denominator. If you put $\rho\approx 10^{-6}$, you're trying to calculate the difference between two numbers that are extremely close to each other, namely, between $1$ and $\exp(\alpha)$ with $\alpha\approx 10^{-21}$. So you run into numerical accuracy problems with machine numbers; see here and at the end of this answer.

In general, when doing numerical calculations for physical problems, it's better to work in dimensionless units, so as to avoid underflow/overflow problems. so if I were you I'd set $h=k_B=c=1$ (this amounts to a choice of units) and go from there.

Answer 2

Intro

Below is a jazzed up version of adaptivePrecisionEvaluate[expr, symbols, precision] I recently used here: Unable to plot $\exp(-a \sqrt{1+x^2})$ for $a>700$.

First, a disclaimer: I figure @acl's advice is better. This merely shows that Mathematica can make the OP's plot. Often advice about numerical scaling is governed by the available numerical tools. I don't expect this will change things, though.

Another reason for the post is the bizarre, randomly different outputs on the same input related to this plot, which is posted at the end.

OP's plot:

bb[T_, h_, rho_, k_, c_] := (2 h rho^3)/c^2 1/(Exp[h rho/(T k)] - 1);

Plot[
 adaptivePrecisionEvaluate[
  bb[5800, 6.626 10^-34, rho, 1.38 10^-23, 299792458]],
 {rho, 0.1 10^-6, 3 10^-6}]

Code for adaptivePrecisionEvaluate[]

ClearAll[adaptivePrecisionEvaluate];
adaptivePrecisionEvaluate // Attributes = {HoldAll};
(*"evaluate a numeric expression:
 *   Uses the adaptive precision functionality of N[expr, prec]
 *   Changes precision of numeric coefficients to infinite
 *   Change down/own/sub/up values of syms to infinite precision
 *   Damn the overhead! Full speed ahead!! ;) "*)
adaptivePrecisionEvaluate[expr_, symbols_ : Automatic, 
   precision_ : Automatic] := With[{
    (* DEFS *)
    approxToExact =(* promote to infinite precision *)
     a : _Real | _Complex :> 
      With[{b = SetPrecision[a, Infinity]}, b /; True],
    syms = Replace[Hold[symbols], {
       Hold[Automatic] :> Thread[
         DeleteDuplicates@
          Cases[Unevaluated[expr], 
           s_Symbol /; 
             MatchQ[Context[Unevaluated@s], "Global`" | Context[]] :> 
            Hold[s], {0, Infinity}, Heads -> True],
         Hold],
       Hold[None] :> Hold[{}],
       Hold[s_List] :> Hold[s],
       Hold[s_Symbol] :> Hold[{s}]
       }],
    prec = Replace[precision, Automatic -> Precision[Hold[expr]]]
    }, With[{
     setExactValues = 
      Function[
       vals, # /. Hold[f_] :> (vals[f] = vals[f] /. approxToExact) &]
     },
    (* MAIN ROUTINE *)
    syms /. Hold[s_] :> Internal`InheritedBlock[s,
       Map[
        (If[System`Private`HasDownEvaluationsQ @@ #,
           setExactValues[DownValues][#],
           If[System`Private`HasOwnEvaluationsQ @@ #,
            setExactValues[OwnValues][#]
            ]
           ];
          If[System`Private`HasSubEvaluationsQ @@ #,
           setExactValues[SubValues][#]
           ];
          If[System`Private`HasUpEvaluationsQ @@ #,
           setExactValues[UpValues][#]
           ]) &
        , Thread[syms]];
       With[{fx = 
          N[Hold[expr] /. approxToExact // ReleaseHold, N@prec]},
        If[prec === MachinePrecision &&
          TrueQ[$MinMachineNumber <= Abs@fx <= $MaxMachineNumber]
         , N[fx]
         , fx]
        ]
       ]]];

Bizarre happenings (V14.0.0, Mac ARM)

If we set the precision to infinity, we get the following version of the OP's function:

adaptivePrecisionEvaluate[
 bb[5800, 6.626  10^-34, rho, 1.38  10^-23, 299792458],
 Automatic, Infinity]
(*
(1065744683248461 rho^3) / 
 (72279169429579204524505290506240000000000000000000000000000000000 *
  (-1 + E^(2486737594246409 rho/300390095145612060000000000000)))
*)

Then the following sometimes yields the same plot above and sometimes a blank plot; it always gives many Power::infy error warnings, no matter which plot is produced. The blank plot is produced more often than the one above, something like 1-10% of the time.

Plot[
 (1065744683248461 rho^3) / 
  (72279169429579204524505290506240000000000000000000000000000000000 *
   (-1 + E^(2486737594246409 rho/300390095145612060000000000000))),
{rho, 0.1  10^-6, 3  10^-6}, WorkingPrecision -> 32]

On WolframCloud Linux, I always get a blank plot.

This evaluation works on both Mac/Linux:

Plot[
 (1065744683248461 rho^3) / 
  (72279169429579204524505290506240000000000000000000000000000000000 *
   (-1 + E^(2486737594246409 rho/300390095145612060000000000000))),
 {rho, N[1/10 10^-6, 32], N[3 10^-6, 32], 1`32*^-7}] // N

I could ask this in a separate Q&A, but since it seems a bug and since it's random, I'm just curious for confirmation:

• Can others replicate this on a Mac? Intel vs. Apple chips?
• What happens on Windows?
• Any reason to think it's not a bug? If so, I'll post a separate Q&A so the answer can get due attention.