2

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 = SystemDumpHeldNumericFunctionQ;

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.

Lauren Pearce
  • 550
  • 2
  • 10
  • 1
    Try replacing expr by Rationalize[Unevaluated@expr, 0] at the end of trappedEvaluation. – Michael E2 Mar 12 '21 at 04:22
  • @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... – Lauren Pearce Mar 12 '21 at 12:42
  • @MichaelE2 If you have a minute to explain why that works, I'd appreciate it. – Lauren Pearce Mar 12 '21 at 12:48
  • @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. – Lauren Pearce Mar 12 '21 at 17:25
  • 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. – Michael E2 Mar 12 '21 at 19:39

0 Answers0