函数插值、求导与积分问题：为何kfm求导结果恒为101.0？

我有三个列表，list_umf是x值，list_kf是y值，list_kfm是list_kf的积分结果（由代码输出得到）。为了验证list_kfm确实是list_kf的积分，我尝试计算list_kfm的导数，理论上结果应该和list_kf一致，但实际计算出的list_kf_re全是101.0，哪里出问题了？

我的代码如下：

import numpy as np
from scipy import integrate, interpolate
from scipy.misc import derivative as deriv
import matplotlib.pyplot as plt

list_kfm = [15.348748494618041, 26.240336614039776, 37.76846357985518, 49.80068952374503, 62.25356792292074, 75.0692188764684, 88.20491343740369, 101.6276911997135,
    115.31128207665246, 129.2342114999071, 143.37856687640036, 157.72915825067278, 172.27292637703843, 186.9985127198004, 201.89593919604192, 216.95636451973587]
list_kf = [168.08871431597626, 179.78615963605742, 188.728883379148, 196.0371678709251, 202.25334207341422, 207.68364358717665, 212.51893919883966, 216.88670040685466,
    220.87653440371076, 224.55397301446894, 227.96847485999652, 231.15833919688876, 234.1538643061246, 236.97945558527186, 239.65507793294745, 242.19728380107006]
list_umf = [0.1, 0.15000000000000002, 0.20000000000000004, 0.25000000000000006, 0.30000000000000004, 0.3500000000000001, 0.40000000000000013, 0.45000000000000007,
    0.5000000000000001, 0.5500000000000002, 0.6000000000000002, 0.6500000000000001, 0.7000000000000002, 0.7500000000000002, 0.8000000000000002, 0.8500000000000002]

f = interpolate.interp1d(
    list_umf, list_kfm, bounds_error=False, fill_value=(15, 217))

list_kf_re = [deriv(f, x) for x in list_umf]

plt.plot(list_umf, list_kfm, label='kfm')
plt.plot(list_umf, list_kf, label='kf')
plt.plot(list_umf, list_kf_re, label='kfre')

print(list_kf_re)
print(list_kf)

问题原因 #

1. 插值与导数步长不匹配：interp1d默认使用线性插值，而scipy.misc.derivative的默认步长dx=1.0过大，计算导数时采样点会超出list_umf的范围，触发你设置的fill_value=(15,217)，最终算出的导数为(217-15)/2=101，完全偏离真实值。
2. 边界填充干扰：bounds_error=False配合不合理的填充值，让插值函数在边界外返回固定值，直接导致导数计算异常。

修复方案 #

方案1：优化插值与导数计算参数 #

使用更高阶的插值方法保证曲线光滑，同时开启边界报错避免填充值干扰，调整导数步长为极小值：

import numpy as np
from scipy import interpolate
from scipy.misc import derivative as deriv
import matplotlib.pyplot as plt

list_kfm = [15.348748494618041, 26.240336614039776, 37.76846357985518, 49.80068952374503, 62.25356792292074, 75.0692188764684, 88.20491343740369, 101.6276911997135,
    115.31128207665246, 129.2342114999071, 143.37856687640036, 157.72915825067278, 172.27292637703843, 186.9985127198004, 201.89593919604192, 216.95636451973587]
list_kf = [168.08871431597626, 179.78615963605742, 188.728883379148, 196.0371678709251, 202.25334207341422, 207.68364358717665, 212.51893919883966, 216.88670040685466,
    220.87653440371076, 224.55397301446894, 227.96847485999652, 231.15833919688876, 234.1538643061246, 236.97945558527186, 239.65507793294745, 242.19728380107006]
list_umf = [0.1, 0.15000000000000002, 0.20000000000000004, 0.25000000000000006, 0.30000000000000004, 0.3500000000000001, 0.40000000000000013, 0.45000000000000007,
    0.5000000000000001, 0.5500000000000002, 0.6000000000000002, 0.6500000000000001, 0.7000000000000002, 0.7500000000000002, 0.8000000000000002, 0.8500000000000002]

# 使用三次插值，开启边界报错避免填充值干扰
f = interpolate.interp1d(list_umf, list_kfm, kind='cubic', bounds_error=True)
# 用极小步长计算导数，避免跨出有效区间
list_kf_re = [deriv(f, x, dx=1e-6) for x in list_umf]

plt.plot(list_umf, list_kfm, label='kfm')
plt.plot(list_umf, list_kf, label='kf')
plt.plot(list_umf, list_kf_re, label='kfre(插值求导)')
plt.legend()
plt.show()

print("插值求导结果：", list_kf_re)
print("原始kf：", list_kf)

方案2：直接使用数值差分（更适合离散数据） #

因为list_kfm是离散积分结果，直接计算相邻点的斜率更准确，无需额外插值：

import numpy as np
import matplotlib.pyplot as plt

list_kfm = [15.348748494618041, 26.240336614039776, 37.76846357985518, 49.80068952374503, 62.25356792292074, 75.0692188764684, 88.20491343740369, 101.6276911997135,
    115.31128207665246, 129.2342114999071, 143.37856687640036, 157.72915825067278, 172.27292637703843, 186.9985127198004, 201.89593919604192, 216.95636451973587]
list_kf = [168.08871431597626, 179.78615963605742, 188.728883379148, 196.0371678709251, 202.25334207341422, 207.68364358717665, 212.51893919883966, 216.88670040685466,
    220.87653440371076, 224.55397301446894, 227.96847485999652, 231.15833919688876, 234.1538643061246, 236.97945558527186, 239.65507793294745, 242.19728380107006]
list_umf = [0.1, 0.15000000000000002, 0.20000000000000004, 0.25000000000000006, 0.30000000000000004, 0.3500000000000001, 0.40000000000000013, 0.45000000000000007,
    0.5000000000000001, 0.5500000000000002, 0.6000000000000002, 0.6500000000000001, 0.7000000000000002, 0.7500000000000002, 0.8000000000000002, 0.8500000000000002]

# 计算相邻点的斜率，补全最后一个点的导数（沿用前一个斜率）
list_kf_diff = np.diff(list_kfm) / np.diff(list_umf)
list_kf_diff = np.append(list_kf_diff, list_kf_diff[-1])

plt.plot(list_umf, list_kfm, label='kfm')
plt.plot(list_umf, list_kf, label='kf')
plt.plot(list_umf, list_kf_diff, label='kfre(数值差分)', linestyle='--')
plt.legend()
plt.show()

print("数值差分结果：", list_kf_diff)
print("原始kf：", list_kf)

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