计算 1/tanh(x) - 1/x 对于非常小的 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

使用 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 2

对于非常小的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 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)

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

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

python

math

numpy

rounding