如何避免 Python 3.7 中 f(x) = (1-cos(x))/x**2 中小数的灾难性取消？

Question

如何避免在 f(x) = (1-cos(x))/x**2 在 Python 3.7?

这是我目前为止尝试过的方法（我知道，关键是一些三角恒等式可以让你避免取消，而且我也知道，使用 L'Hopital 规则，f 的极限→0 (x) 是 0.5，所以正确的程序输出是非常接近 0.5 的东西，这就是你使用 x = 1.2e-4 时得到的结果，但是你会用较小的数字（如 1.2e-8）取消，并且我需要做到这一点，这样才不会发生。

from math import *
def f(x):     #these are all the same function using different identities  
   a = (1-(sin(x)/tan(x)))/(x**2)
   b = (1-(sin(2*x)/(2*sin(x))))/(x**2)
   c = (1-((1-((tan(x/2))**2))/(1+(tan(x/2))**2)))/(x**2)
   d = (sin(x)**2+cos(x)**2-cos(x))/(x**2)
   e = (sin(x)**2+cos(x)**2-(sin(2*x)/(2*sin(x))))/(x**2)
   return a, b, c, d, e

print(k(1.2e-8))
#Output: (0.0, 0.7709882115452477, 0.0, 0.0, 0.0) - whereas need 0.5000...

Answer 1

像这样：

sin(x)/x * tan(x/2)/x

一直到最后，x = 1e-308 还可以。

不幸的是，我无法深入了解为什么它运作良好。

Answer 2

问题是 float 和 double 的精度有限。您必须使用更精确的算术，例如 mpfr。可以通过绑定如

在Python中使用

https://pypi.org/project/gmpy2/

这是一个示例，我在一个名为 Sagemath 的更高级环境中使用它：我使用的是 100 位精度：

sage: R = RealField(100)  // 100 bits of precision
sage: def f(x):     #these are all the same function using different identities  
....:    a = (1-(sin(x)/tan(x)))/(x**2)
....:    b = (1-(sin(2*x)/(2*sin(x))))/(x**2)
....:    c = (1-((1-((tan(x/2))**2))/(1+(tan(x/2))**2)))/(x**2)
....:    d = (sin(x)**2+cos(x)**2-cos(x))/(x**2)
....:    e = (sin(x)**2+cos(x)**2-(sin(2*x)/(2*sin(x))))/(x**2)
....:    return a, b, c, d, e
....: 
sage: f(R(1.2e-8))
(0.50000000000000264647624775223,
 0.49999999999999716827551705076,
 0.49999999999999716827551705076,
 0.49999999999999716827551705076,
 0.49999999999999169007478634929)

Answer 3

改用numpy.float128。 numpy 是数据分析和更复杂数学的标准。它可以在终端中使用以下命令安装。

pip install numpy

from numpy import *
def f(x):     #these are all the same function using different identities  
   a = (1-(sin(x)/tan(x)))/(x**2)
   b = (1-(sin(2*x)/(2*sin(x))))/(x**2)
   c = (1-((1-((tan(x/2))**2))/(1+(tan(x/2))**2)))/(x**2)
   d = (sin(x)**2+cos(x)**2-cos(x))/(x**2)
   e = (sin(x)**2+cos(x)**2-(sin(2*x)/(2*sin(x))))/(x**2)
   return a, b, c, d, e
print(f(float128(1.2e-8)))

这会打印

(0.5003141275115400736, 0.49956120933620291774, 0.49993766842387149567, 0.49993766842387149567, 0.49956120933620291774)

Answer 4

您可以有条件地 return 小 x 的限制，如下所示：

import matplotlib.pyplot as plt
from numpy import *
epsilon=1e-8

def f(x):
   if x<epsilon
      return 0.5

   return (1-sin(x)/tan(x))/x**2
   #note: same as (1-cos(x))/x**2

x=arange(0,6,0.01)
y=vectorize(f)
plt.plot(x,y(x))
plt.show()

绘制的曲线看起来很平滑

注意：比起数学，我更喜欢 numpy。 vectorize使得用数组调用函数成为可能（效率不是很高，但是好用）。

如何避免 Python 3.7 中 f(x) = (1-cos(x))/x**2 中小数的灾难性取消？

How do I avoid catastrophic cancellation for small numbers in f(x) = (1-cos(x))/x**2 in Python 3.7?

python

trigonometry

function

floating-accuracy

python-3.x