1

How do I evaluate xe^x/(e^x-1) with numerical stability around zero and when x is very positive or negative? I have access to all the usual mathematical functions in numpy and scipy.

Neil G
  • 32,138
  • 39
  • 156
  • 257
  • Just found scipy.special.exprel, which would work. I need to have my code work on theano or tensorflow as well, so a more basic solution would be better. – Neil G Sep 03 '16 at 04:03
  • Frustratingly, `np.expm1` seems to give `nan` instead of `inf` for large inputs. That makes it harder to avoid things like `if` or `where`. – user2357112 Sep 03 '16 at 04:33
  • This seems to be a [platform-dependent bug](https://github.com/numpy/numpy/issues/6818), still open. – user2357112 Sep 03 '16 at 04:35

1 Answers1

2
def f(x):
  if abs(x) > 0.1: return x*exp(x)/(exp(x)-1)
  else: return 1/(1.-x/2.+x**2/6.-x**3/24.)

The expansion in the last line can be extended in the obvious fashion if more precision is required, and can be made faster by re-phrasing. As it stands, it errs by as much as 1e-6.

P. Pearson
  • 116
  • 2
  • 2
    There is a not so large value `z` such that for `x > z` the real value of the expression `x*exp(x)/(exp(x)-1)` is very close to `x`, whereas computing it straightforwardly will result in overflow and nan. – Leon Sep 03 '16 at 05:49
  • Good point, @Leon. To handle this possibility, the cases of `x > 0.1` and `x < -0.1` should be handled separately. If `x > 0.1`, return `x/(1-exp(-x))`; if `x < -0.1`, return `x*exp(x)/(exp(x)-1)`; and otherwise use the polynomial approximation. – P. Pearson Sep 08 '16 at 03:49