使用SymPy计算20阶以上高阶导数结果异常问题咨询

问题背景 #

我正在利用Sellmeier方程研究介质的高阶色散，定义了以下两个表达式：

n_expr = sp.sqrt( 7+ 8/λ + 9/(λ**2 - 5) + 5/(6-λ**2) )
λ_ω_expr = c/ω

需要计算复合函数 k_expr = ω * n_expr.subs(λ, λ_ω_expr) / c 对ω的1至30阶导数。但直接使用sp.diff(k_expr, ω, n)求值的结果，与手动通过链式法则分步计算的结果不一致，且根据研究背景判断，手动链式法则的结果是正确的。手动定义30阶导数过于繁琐，因此需要定位直接求导方法出错的原因。

最小可复现代码 #

通用代码 #

import sympy as sp

f = (75 - 24.5) * (75 + 570.83)
c = 2.99792458e11

λ = sp.Symbol('λ')
ω = sp.Symbol('ω')

ω_test = 1.2e15

n_expr = sp.sqrt(
    5.756
    + 2.86e-6 * f
    + (0.0983 + 4.7e-8 * f) / (λ**2 - (0.202 + 6.113e-8 * f)**2)
    + (189.32 + 1.516e-4 * f) / (λ**2 - 12.52**2)
    - 1.32e-2 * λ**2
)

n_func = sp.lambdify(λ, n_expr, "numpy", cse=True)

λ_ω_expr = c/ω

手动链式法则（结果正确） #

dλ_dω_expr      = sp.diff(λ_ω_expr, ω)
d2λ_dω2_expr    = sp.diff(λ_ω_expr, ω, 2)
dn_dλ_expr = sp.diff(n_expr, λ)
d2n_dλ2_expr = sp.diff(n_expr, λ, 2)

dλ_dω_func      = sp.lambdify(ω, dλ_dω_expr, "numpy")
d2λ_dω2_func    = sp.lambdify(ω, d2λ_dω2_expr, "numpy")
dn_dλ_func = sp.lambdify(λ, dn_dλ_expr, "numpy")
d2n_dλ2_func = sp.lambdify(λ, d2n_dλ2_expr, "numpy")

def correct_dk_dω(ω):
    return ( 1*n_func(c/ω) + ω*dn_dλ_func(c/ω)*dλ_dω_func(ω) )/c

def correct_d2k_dω2(ω):
    return ( 2*dn_dλ_func(c/ω)*dλ_dω_func(ω) + ω*d2n_dλ2_func(c/ω)*d2λ_dω2_func(ω) )/c

print('correct 1order:', correct_dk_dω(ω_test))
print('correct 2order:', correct_d2k_dω2(ω_test))

直接求导方法（结果错误） #

k_expr = ω * n_expr.subs(λ, λ_ω_expr) / c

wrong_dk_dω_expr = sp.diff(k_expr, ω)
wrong_d2k_dω2_expr = sp.diff(k_expr, ω, 2)
wrong_dk_dω_func = sp.lambdify(ω, wrong_dk_dω_expr, "numpy")
wrong_d2k_dω2_func = sp.lambdify(ω, wrong_d2k_dω2_expr, "numpy")

print('wrong 1order:', wrong_dk_dω_func(ω_test))
print('wrong 2order:', wrong_d2k_dω2_func(ω_test))

可能的原因分析 #

• 符号替换与化简顺序问题：直接替换λ = c/ω后再求导，SymPy的自动化简可能引入了符号层面等价但数值计算时精度丢失的表达式，或者默认了某些分母非零的假设，与实际求值场景冲突。
• 复合函数求导的展开差异：手动链式法则是分步求导后再代入数值，而直接求导是先合成复杂表达式再求导，生成的表达式包含大量冗余项，在转换为数值函数时，浮点数运算的累积误差被放大，尤其是高阶导数对数值误差极为敏感。
• SymPy化简策略的影响：直接求导后的表达式未完全展开或合并同类项，lambdify转换时生成了低效的数值运算代码，进一步加剧精度损失。可尝试在求导后调用sp.simplify()或sp.expand()再转换。
• lambdify的转换精度问题：复杂表达式转换为numpy函数时，部分符号运算的转换逻辑可能存在精度损耗，而手动分步转换的函数结构更简单，数值计算稳定性更高。

验证建议 #

1. 对直接求导后的表达式执行符号化简：

   simplified_dk_dω = sp.simplify(wrong_dk_dω_expr)
   simplified_func = sp.lambdify(ω, simplified_dk_dω, "numpy")
   

   再对比求值结果。
2. 检查两种方法的表达式是否符号等价：

   # 先手动推导正确的符号表达式，再做差化简
   correct_dk_dω_expr_symbolic = (n_expr.subs(λ, c/ω) + ω * dn_dλ_expr.subs(λ, c/ω) * dλ_dω_expr)/c
   diff_expr = sp.simplify(wrong_dk_dω_expr - correct_dk_dω_expr_symbolic)
   

   若结果为0，说明是数值精度问题；若不为0，说明SymPy求导过程存在错误。

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