Python实现导数计算器遇h<1e-7时精度丢失问题求助

• 浮点数舍入误差：Python中的float属于IEEE 754双精度浮点数，仅能保留约15-17位有效数字。当h过小（如小于1e-7）时，x+h与x的浮点数表示会因精度限制发生舍入，导致f(x+h)-f(x)的计算出现严重误差。以x=7为例，双精度浮点数在7附近的最小可分辨增量（ULP，单位最后位置）约为1.55e-15，当h=1e-14时，7+1e-14的浮点数结果会被舍入为近似值，而非精确的7+1e-14，这使得差值计算的误差被放大。
• 误差平衡被打破：向前差分公式(f(x+h)-f(x))/h的截断误差随h增大而增大（O(h)量级），舍入误差随h减小而增大。当h小于临界值（约1e-7）时，舍入误差的增长幅度超过截断误差的减少幅度，整体精度开始骤降。

1. 替换为中心差分公式 #

中心差分的截断误差为O(h²)，精度远高于向前差分，同时能有效缓解舍入误差的影响。公式为：

def central_diff(x, h):
    return (f(x + h) - f(x - h)) / (2 * h)

对于f(x)=x²，中心差分的计算结果可简化为2x，理论上无截断误差，实际计算中受浮点数精度影响极小。

2. 选择最优h值 #

对于双精度浮点数，最优h通常取√(ε)（ε为机器精度，Python中可通过float.epsilon获取，约2.22e-16），即h≈1e-8。这个值能平衡截断误差与舍入误差，使整体计算精度达到最高。

3. 优化后的完整代码示例 #

from math import sqrt

def f(x):
    return x ** 2

def central_diff(x, h=None):
    if h is None:
        # 自动计算最优h值，适配不同大小的x
        h = sqrt(float.epsilon) * max(1.0, abs(x))
    return (f(x + h) - f(x - h)) / (2 * h)

a = float(input("At what point? "))
print("高精度导数结果:", central_diff(a))

# 可选：对比不同h值的计算效果
for exp in range(-1, -16, -1):
    h = 10 ** exp
    print(f"h=1e{exp}: {central_diff(a, h)}")

4. 进阶方案：自动微分 #

若需彻底避免数值误差，可使用自动微分方法——通过跟踪计算过程中的导数传递，直接得到解析导数结果。示例如下：

from autograd import grad

def f(x):
    return x ** 2

# 生成f(x)的导数函数
df = grad(f)
print(df(7.0))  # 精确输出14.0

内容的提问来源于stack exchange，提问作者theauthor