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.
exprbyRationalize[Unevaluated@expr, 0]at the end oftrappedEvaluation. – Michael E2 Mar 12 '21 at 04:22trappedEvaluation[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:42Plot[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:25trappedEvaluation[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], becauseexprwas evaluating to machine precision. I thought adding an extraRationalizemight fix it. The problem withPlotseems to be that this leads toRationalize[f[x]]with a symbolicxinstead of a number. Afterwardsxevaluates to a machine-precision number. Another potential issue is thatPlotcan plot only machine real coordinates. So between around 10^±308 in magnitude. – Michael E2 Mar 12 '21 at 19:39