我一直在Numpy/Scipy中寻找包含有限差分函数的模块。但是，我找到的最接近的东西是numpy.gradient()，它对于2阶精度的1阶有限差分很有用，但是如果您想要更高阶的导数或更精确的方法，那么就没有那么多了。对于这种事情，我什至还没有找到很多特定的模块。大多数人在需要时似乎正在做“自己动手”的事情。所以我的问题是，是否有人知道专门针对高阶(精度和微分)有限差分方法的任何模块(Numpy/Scipy的一部分或第三方模块)。我有自己的代码正在使用，但是目前有点慢，如果已有可用的东西，我不会尝试对其进行优化。
请注意，我在说的是有限差分，而不是导数。我已经看过scipy.misc.derivative()和Numdifftools，它们都是我没有的解析函数的派生类。
最佳答案
快速做到这一点的一种方法是通过与高斯核的导数进行卷积。最简单的情况是数组与[-1, 1]的卷积，给出了简单的有限差分公式。除此之外，(f*g)'= f'*g = f*g'是卷积的*，因此您最终将派生函数与普通高斯卷积在一起，因此当然可以稍微平滑数据，可以通过选择最小的合理内核来使其最小化。
import numpy as np
from scipy import ndimage
import matplotlib.pyplot as plt
#Data:
x = np.linspace(0,2*np.pi,100)
f = np.sin(x) + .02*(np.random.rand(100)-.5)
#Normalization:
dx = x[1] - x[0] # use np.diff(x) if x is not uniform
dxdx = dx**2
#First derivatives:
df = np.diff(f) / dx
cf = np.convolve(f, [1,-1]) / dx
gf = ndimage.gaussian_filter1d(f, sigma=1, order=1, mode='wrap') / dx
#Second derivatives:
ddf = np.diff(f, 2) / dxdx
ccf = np.convolve(f, [1, -2, 1]) / dxdx
ggf = ndimage.gaussian_filter1d(f, sigma=1, order=2, mode='wrap') / dxdx
#Plotting:
plt.figure()
plt.plot(x, f, 'k', lw=2, label='original')
plt.plot(x[:-1], df, 'r.', label='np.diff, 1')
plt.plot(x, cf[:-1], 'r--', label='np.convolve, [1,-1]')
plt.plot(x, gf, 'r', label='gaussian, 1')
plt.plot(x[:-2], ddf, 'g.', label='np.diff, 2')
plt.plot(x, ccf[:-2], 'g--', label='np.convolve, [1,-2,1]')
plt.plot(x, ggf, 'g', label='gaussian, 2')
由于您提到
np.gradient，所以我假设您至少有2d数组，因此以下内容适用于此:如果要对ndarrays这样做，它内置在scipy.ndimage包中。但是请谨慎，因为这当然不能给您完整的渐变，但是我相信各个方向的乘积。希望拥有更好专业知识的人会说出来。这是一个例子:
from scipy import ndimage
x = np.linspace(0,2*np.pi,100)
sine = np.sin(x)
im = sine * sine[...,None]
d1 = ndimage.gaussian_filter(im, sigma=5, order=1, mode='wrap')
d2 = ndimage.gaussian_filter(im, sigma=5, order=2, mode='wrap')
plt.figure()
plt.subplot(131)
plt.imshow(im)
plt.title('original')
plt.subplot(132)
plt.imshow(d1)
plt.title('first derivative')
plt.subplot(133)
plt.imshow(d2)
plt.title('second derivative')
使用
gaussian_filter1d可使您沿特定轴取方向导数:imx = im * x
d2_0 = ndimage.gaussian_filter1d(imx, axis=0, sigma=5, order=2, mode='wrap')
d2_1 = ndimage.gaussian_filter1d(imx, axis=1, sigma=5, order=2, mode='wrap')
plt.figure()
plt.subplot(131)
plt.imshow(imx)
plt.title('original')
plt.subplot(132)
plt.imshow(d2_0)
plt.title('derivative along axis 0')
plt.subplot(133)
plt.imshow(d2_1)
plt.title('along axis 1')
上面的第一组结果令我有些困惑(当曲率指向下方时，原始峰显示为二阶导数的峰值)。在不进一步研究2d版本如何工作的情况下，我只能真正推荐1d版本。如果要幅度，只需执行以下操作:
d2_mag = np.sqrt(d2_0**2 + d2_1**2)