Python:任意精度

Question

我尝试制作这个简单的程序，以正常精度计算函数的导数：

# second derivative of a function
def diff2(f, x, h=1E-6):
        r = (f(x-h) - 2*f(x) + f(x+h))/float(h*h)
        return r
# define the function to derivate
def g(t):
        return t**(-6)
# decresing h increasing the precision of the derivative
# ROUND-OFF problems are present
for k in range(1,15):
        h = 10**(-k) # 15 different value of h
        d2g = diff2(g, 1, h) # compute d'' of g 15-th times in point t=1
        print 'h=%.0e: %.5f' % (h, d2g)

从打印操作中可以看出，当 k 由于四舍五入而大于 8 时，我遇到了问题。我知道我可以使用：

from decimal import *

但是我不知道如何在我的函数中实现这些命令。

有人可以帮我吗？

Answer 1

值得研究一下python模块mpmath，它可以处理任意精度。例如：

>>> from mpmath import mp
>>> mp.dps = 50
>>> print(mp.quad(lambda x: mp.exp(-x**2), [-mp.inf, mp.inf]) ** 2)
3.1415926535897932384626433832795028841971693993751

您可以简单地更改类型，让您的函数更精确地工作。然而，值得注意的是@halex的评论和答案。

Answer 2

如果您希望导

数位于位置x_0并使用浮点数执行计算，则最小化数值误差的h的最佳值是sqrt(sys.float_info.epsilon)*x_0对于x_0=1的情况，该值约为 1E-8。

有关此值的更多信息和推导，请参阅第 How to Choose h 章，从第 4 页开始，直到这个关于数值微分的简短脚本的结尾。

Answer 3

您可以使用十进制模块：

from decimal import Decimal
# second derivative of a function
def diff2(f, x, h=1E-6):
    x, h = Decimal(x), Decimal(h)
    r = (f(x - h) - 2 * f(x) + f(x + h)) / Decimal(h * h)
    return r