Evaluate 1/tanh(x) - 1/x for very small x

Question

I need to compute the quantity

1/tanh(x) - 1/x

for x > 0, where x can be both very small and very large.

Asymptotically for small x, we have

1/tanh(x) - 1/x  ->  x / 3

and for large x

1/tanh(x) - 1/x  ->  1

Anyhow, when computing the expression, already from 10^-7 and smaller round-off errors lead to the expression being evaluated as exactly 0:

import numpy
import matplotlib.pyplot as plt


x = numpy.array([2**k for k in range(-30, 30)])
y = 1.0 / numpy.tanh(x) - 1.0 / x

plt.loglog(x, y)
plt.show()

Answer 1

For very small x, one could use the Taylor expansion of 1/tanh(x) - 1/x around 0,

y = x/3.0 - x**3 / 45.0 + 2.0/945.0 * x**5

The error is of the order O(x**7), so if 10^-5 is chosen as the breaking point, relative and absolute error will be well below machine precision.

import numpy
import matplotlib.pyplot as plt


x = numpy.array([2**k for k in range(-50, 30)])

y0 = 1.0 / numpy.tanh(x) - 1.0 / x
y1 = x/3.0 - x**3 / 45.0 + 2.0/945.0 * x**5
y = numpy.where(x > 1.0e-5, y0, y1)


plt.loglog(x, y)
plt.show()

Answer 2

A probably simpler solution to overcome this is changing the data type under which numpy is operating:

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(-30, 30, dtype=np.longdouble)
x = 2**x
y = 1.0 / np.tanh(x) - 1.0 / x

plt.loglog(x, y)
plt.show()

Using longdouble as data type does give the proper solution without rounding errors.

---

I did sightly modify your example, in your case the only thing you need to modify is:

x = numpy.array([2**k for k in range(-30, 30)])

to:

x = numpy.array([2**k for k in range(-30, 30)], dtype=numpy.longdouble)

Answer 3

Use the python package mpmath for arbitrary decimal precision. For example:

import mpmath
from mpmath import mpf

mpmath.mp.dps = 100 # set decimal precision

x = mpf('1e-20')

print (mpf('1') / mpmath.tanh(x)) - (mpf('1') / x)
>>> 0.000000000000000000003333333333333333333333333333333333333333311111111111111111111946629156220629025294373160489201095913

It gets extremely precise.

Look into mpmath plotting. mpmath plays well with matplotlib, which you are using, so this should solve your problem.

---

Here is an example of how to integrate mpmath into the code you wrote above:

import numpy
import matplotlib.pyplot as plt
import mpmath
from mpmath import mpf

mpmath.mp.dps = 100 # set decimal precision

x = numpy.array([mpf('2')**k for k in range(-30, 30)])
y = mpf('1.0') / numpy.array([mpmath.tanh(e) for e in x]) - mpf('1.0') / x

plt.loglog(x, y)
plt.show()