2
$\begingroup$

In contrast to earlier versions of Mathematica, in versions 11.3 and later Mathematica no longer switches to arbitrary precision when it detects a MachineNumber underflow. While this behavior is consistent with IEEE 754, this can also be problematic.

Background (and a workaround for the exponential function can be found here: New General::munfl error and loss of precision However, since functions other than exponentials can cause underflows, a more general workaround is useful. I am trying to use the method proposed in the second solution of: Underflow error General::munfl from E^x instead of Exp[x]

However, the trappedEvaluation function defined there is failing to detect or fix the underflow in the following:

(* This is the function from the above-mentioned solution *)
SetAttributes[trappedEvaluation, HoldAll]

trappedEvaluation[expr_] := Check[
   expr, 
   numericize@Block[{Real}, 
       f_?nfQ[a___, r_Real, b___] ^:= Rationalize[Unevaluated@f[a, r, b], 0];
       expr
   ], 
   General::munfl
 ]

 nfQ = System`Dump`HeldNumericFunctionQ;

 numericize[e_] :=  With[{res = N[e, 16]}, 
     Quiet@Check[
        N[res], 
        res, 
        General::munfl
     ]
  ]

 (* Here is my function: *)
 f[x_] = Exp[-(1/x)]*WhittakerW[-(I/x), 1/4, I*x];

 (* Test: *)
 trappedEvaluation[f[.0001]]

This still returns 0.+0.i. In fact, based on the analytic expansions in the DLMF this should grow exponentially as x approaches 0, and the behavior up until the underflow error is consistent with this.

Ultimately, I need to make a log plot of f[x]*Conjugate[f[x]], and f'[x]*Conjugate[f'[x]] so the solution has be able to work at the arbitrary points that Mathematica chooses when plotting. Therefore rationalizing by hand (or setting the precision by hand) isn't a solution. I'm afraid the derivative will also cause issues, but first I want to get it working on my function.

I don't believe the problem is in my implementation of trappedEvaluation, as it works appropriately on the function in Underflow error General::munfl from E^x instead of Exp[x]. Namely,

   sigma = 0.01;
   gaussian = Function[x, 1/(Sqrt[2 Pi] sigma) Exp[-1/2 (x/sigma)^2]]
   gaussianP[x_] = D[gaussian[x], x]
   trappedEvaluation[gaussianP[1.5]]

returns -9.206189566469410*10^-4881 as it should. I would much appreciate it if someone could help me understand how to generalizing the trappedEvaluation function appropriately.

$\endgroup$
5
  • 1
    $\begingroup$ Try replacing expr by Rationalize[Unevaluated@expr, 0] at the end of trappedEvaluation. $\endgroup$ – Michael E2 Mar 12 at 4:22
  • $\begingroup$ @MichaelE2 I realized you probably meant to change: trappedEvaluation[expr_] := Check[expr, numericize@ Block[{Real}, f_?nfQ[a___, r_Real, b___] ^:= Rationalize[Unevaluated@f[a, r, b], 0]; expr/. {expr -> Rationalize[Unevaluated@expr, 0]}], General::munfl] so I made that change instead. That seems to work... $\endgroup$ – Lauren Pearce Mar 12 at 12:42
  • $\begingroup$ @MichaelE2 If you have a minute to explain why that works, I'd appreciate it. $\endgroup$ – Lauren Pearce Mar 12 at 12:48
  • $\begingroup$ @MichaelE2 : PS- it appears not to work if called within plot (e.g., Plot[Re[trappedEvaluation[f[x]]],{x,0,0.001}] appears to plot zero uniformly. Obviously an interpolated function can be built, but it'd be nice if I could understand why what you did work to see if I can modify it. If you could explain I'd much appreciate it. $\endgroup$ – Lauren Pearce Mar 12 at 17:25
  • $\begingroup$ I had in mind trappedEvaluation[expr_] := Check[expr, numericize@Block[{Real}, f_?nfQ[a___, r_Real, b___] ^:= Rationalize[Unevaluated@f[a, r, b], 0]; Rationalize[Unevaluated@expr, 0]], General::munfl], because expr was evaluating to machine precision. I thought adding an extra Rationalize might fix it. The problem with Plot seems to be that this leads to Rationalize[f[x]] with a symbolic x instead of a number. Afterwards x evaluates to a machine-precision number. Another potential issue is that Plot can plot only machine real coordinates. So between around 10^±308 in magnitude. $\endgroup$ – Michael E2 Mar 12 at 19:39

Your Answer

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

Browse other questions tagged or ask your own question.