为何Scipy计算的函数导数随dx参数变化？如何验证结果正确性？

问题背景 #

尝试用SymPy和Scipy计算同一函数的四阶导数，结果差异悬殊：

• SymPy计算结果：-73035.8044625845
• Scipy（dx=1e-6）计算结果：15154544286.133389

使用的代码如下：

import sympy as smp
import numpy as np
import scipy as sp

# SymPy符号求导
x = smp.symbols('x', real=True)
f = smp.exp(-smp.sin(x**2)) * smp.sin(2**x) * smp.log(3*smp.sin(x)**2/x)
smpval = smp.diff(f, x, 4).subs([(x, 4)]).evalf()

# Scipy数值求导
def f(x):
    return sp.exp(-sp.sin(x**2)) * sp.sin(2**x) * sp.log(3*sp.sin(x)**2/x)

from scipy.misc import derivative
spval = derivative(f, x0=4, dx=1e-6, n=4, order=5)
print(smpval)
print(spval)

偏差原因 #

1. 计算原理本质不同 #

• SymPy采用符号计算：先对函数进行精确的符号推导，得到四阶导数的解析表达式后，再代入x=4做浮点求值，全程仅在最后一步存在浮点精度误差，结果准确可靠。
• Scipy的derivative采用数值差分近似：基于泰勒展开的有限差分公式近似导数，属于近似计算，高阶导数对参数极其敏感。

2. 高阶数值导数的固有不稳定性 #

四阶导数的5点差分公式为：

f''''(x0) ≈ [f(x0-2dx) - 4f(x0-dx) + 6f(x0) -4f(x0+dx) + f(x0+2dx)] / (dx⁴)

当dx=1e-6时，分母为1e-24，分子中微小的浮点舍入误差会被放大1e24倍，直接导致结果完全失真，这就是你得到异常大数值的核心原因。

3. 函数特性加剧误差 #

你的函数包含sin(2^x)（高频振荡）、log(3sin²x/x)（局部变化剧烈）等成分，在x=4附近函数变化率极高，进一步放大了数值差分的误差。

验证方法 #

方法一：以SymPy结果为基准 #

SymPy的符号求导是精确推导，只要函数定义无误，结果就是浮点精度内的准确值。可以将SymPy导出的四阶导数表达式转换为数值函数，代入x=4计算，与原结果对比验证：

# 将SymPy的四阶导数转为数值函数
f4_sym = smp.lambdify(x, smp.diff(f, x, 4), 'numpy')
print(f4_sym(4))  # 应与smpval一致

方法二：优化Scipy计算参数 #

调整dx到合理范围（如1e-3），同时增大差分模板的阶数order（如9），让数值结果收敛到SymPy的数值：

spval_optimized = derivative(f, x0=4, dx=1e-3, n=4, order=9)
print(spval_optimized)  # 结果会接近-73000左右

方法三：分步验证低阶导数 #

先对比一阶、二阶导数的SymPy与Scipy结果（调整dx到合适值），确认低阶时两者一致，再推导高阶导数的正确性，间接验证四阶导数结果。

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