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
    Commented Mar 12, 2021 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$ Commented Mar 12, 2021 at 12:42
  • $\begingroup$ @MichaelE2 If you have a minute to explain why that works, I'd appreciate it. $\endgroup$ Commented Mar 12, 2021 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$ Commented Mar 12, 2021 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
    Commented Mar 12, 2021 at 19:39

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.