对非常小的 x 计算 1/tanh(x) - 1/x

Question

我需要计算数量

1/tanh(x) - 1/x

对于x > 0，其中x 可以非常小也可以非常大。

对于小的x，我们有

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

对于大型x

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

无论如何，在计算表达式时，已经从 10^-7 和较小的舍入错误导致表达式被精确地评估为 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

对于非常小的x，可以使用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

误差是O(x**7)的顺序，所以如果选择10^-5作为断点，相对误差和绝对误差将远低于机器精度。

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

使用 python 包mpmath 获得任意小数精度。例如：

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

它变得非常精确。

查看mpmath plotting。 mpmath 与您正在使用的 matplotlib 配合得很好，所以这应该可以解决您的问题。

---

这是一个如何将 mpmath 集成到您上面编写的代码中的示例：

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()

Answer 3

一个可能更简单的解决方案是改变 numpy 运行的数据类型：

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()

使用longdouble 作为数据类型确实可以提供正确的解决方案而不会出现舍入错误。

---

我确实修改了您的示例，在您的情况下，您唯一需要修改的是：

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

到：

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