Python 有限差分函数？

Question

我一直在 Numpy/Scipy 中寻找包含有限差分函数的模块。 然而，我发现的最接近的是

numpy.gradient()

，它对于二阶精度的一阶有限差分很有用，但如果您想要更高阶的导数或更准确的方法，则效果就不那么好了。 我什至还没有找到很多专门用于此类事情的模块；大多数人似乎都在根据需要做“自己动手”的事情。 所以我的问题是，是否有人知道专门用于高阶（精度和导数）有限差分方法的模块（Numpy/Scipy 的一部分或第三方模块）。 我正在编写自己的代码，但目前速度有点慢，如果已经有可用的东西，我不会尝试优化它。

请注意，我说的是有限差分，而不是导数。 我见过

scipy.misc.derivative()

和 Numdifftools，它们采用分析函数的导数，但我没有。

Answer 1

快速做到这一点的一种方法是与高斯核的导数进行卷积。 简单的情况是数组与

[-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 数组，所以以下内容适用于此：如果你想对 ndarray 执行此操作，则将其内置于

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)

Answer 2

绝对喜欢skewchan给出的答案。 这是一项很棒的技术。 但是，如果您需要使用

numpy.convolve

我想提供这个小修复。 而不是做：

#First derivatives:
cf = np.convolve(f, [1,-1]) / dx
....
#Second derivatives:
ccf = np.convolve(f, [1, -2, 1]) / dxdx
...
plt.plot(x, cf[:-1], 'r--', label='np.convolve, [1,-1]')
plt.plot(x, ccf[:-2], 'g--', label='np.convolve, [1,-2,1]')

...使用

'same'

中的

numpy.convolve

选项，如下所示：

#First derivatives:
cf = np.convolve(f, [1,-1],'same') / dx
...
#Second derivatives:
ccf = np.convolve(f, [1, -2, 1],'same') / dxdx
...
plt.plot(x, cf, 'rx', label='np.convolve, [1,-1]')
plt.plot(x, ccf, 'gx', label='np.convolve, [1,-2,1]')

...以避免差一索引错误。

绘图时还要注意 x 索引。

numy.diff

和

numpy.convolve

的点必须相同！ 要修复差一错误（使用我的

'same'

代码），请使用：

plt.plot(x, f, 'k', lw=2, label='original')
plt.plot(x[1:], df, 'r.', label='np.diff, 1')
plt.plot(x, cf, 'rx', label='np.convolve, [1,-1]')
plt.plot(x, gf, 'r', label='gaussian, 1')
plt.plot(x[1:-1], ddf, 'g.', label='np.diff, 2')
plt.plot(x, ccf, 'gx', label='np.convolve, [1,-2,1]')
plt.plot(x, ggf, 'g', label='gaussian, 2')

您可能想看看 findiff 项目。我自己尝试过，它可以让您方便地对任何维度、任何导数阶数和任何所需精度阶数的 numpy 数组求导数。该项目网站称其特点：

• 沿任意轴对任意维数的数组进行微分
• 任何所需阶数的偏导数
• 向量微积分中的标准运算符，如梯度、散度和旋度
• 可以处理均匀和非均匀网格
• 可以处理具有常数和可变系数的导数的任意线性组合
• 可指定精度顺序
• 完全矢量化以提高速度
• 计算均匀和非均匀网格的任何阶数和精度的原始有限差分系数

Answer 3

另一种方法是对数据进行微分插值。这是由 unutbu 建议的，但我没有看到此处或任何链接问题中使用的方法。例如，

UnivariateSpline

，来自

scipy.interpolate

有一个有用的内置导数方法。

import numpy as np
from scipy.interpolate import UnivariateSpline
import matplotlib.pyplot as plt

# data
n = 1000
x = np.linspace(0, 100, n)
y = 0.5 * np.cumsum(np.random.randn(n))

k = 5 # 5th degree spline
s = n - np.sqrt(2*n) # smoothing factor
spline_0 = UnivariateSpline(x, y, k=k, s=s)
spline_1 = UnivariateSpline(x, y, k=k, s=s).derivative(n=1)
spline_2 = UnivariateSpline(x, y, k=k, s=s).derivative(n=2)

# plot data, spline fit, and derivatives
fig = plt.figure()
ax = fig.add_subplot(111)

ax.plot(x, y, 'ko', ms=2, label='data')
ax.plot(x, spline_0(x), 'k', label='5th deg spline')
ax.plot(x, spline_1(x), 'r', label='1st order derivative')
ax.plot(x, spline_2(x), 'g', label='2nd order derivative')

ax.legend(loc='best')
ax.grid()

注意样条曲线拟合（黑色曲线）的波峰和波谷处的一阶导数（红色曲线）的零交叉。

Answer 4

FastFD 可以提供帮助：

pip install fastfd

基本文档可以在这里找到：

https://github.com/stefanmeili/FastFD

这里：

https://github.com/stefanmeili/FastFD/tree/main/docs/examples

Answer 5

更新了原始函数和用不同方法计算的导数的图：

import numpy as np
from scipy import ndimage
import matplotlib.pyplot as plt
from scipy.interpolate import splrep, splev

#Data:
n = 100
x = np.linspace(0,2*np.pi,n)
f = np.sin(x) + .02*(np.random.rand(n)-.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
tck = splrep(x, f, k=2)
sf1 = splev(x, tck, der=1)
ng1 = np.gradient(f, 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
sf2 = splev(x, tck, der=2)
ng2 = np.gradient(ng1, dx)

#Plotting:
fig, ax = plt.subplots(layout='constrained', figsize=(15, 8))
ax.plot(x, f, 'k', lw=2, label='original')
ax.set_ylabel('original function')

ax2 = ax.twinx()
ax2.plot(x[:-1], df, 'g', label='np.diff, 1')
ax2.plot(x, cf[:-1], 'r', label='np.convolve, [1,-1]')
ax2.plot(x, gf, 'y', label='gaussian, 1')
ax2.plot(x, sf1, 'b', label='splev der=1')
ax2.plot(x, ng1, 'm', label='numpy.gradient 1')
ax2.set_ylabel('first derivative')

plt.title('Original data and 1st derivatives')
plt.legend()
plt.grid()


#Plotting:
fig, ax = plt.subplots(layout='constrained', figsize=(15, 8))
ax.plot(x, f, 'k', lw=2, label='original')
ax.set_ylabel('original function')

ax2 = ax.twinx()
ax2.plot(x[:-2], ddf, 'g', label='np.diff, 2')
ax2.plot(x, ccf[:-2], 'r', label='np.convolve, [1,-2,1]')
ax2.plot(x, ggf, 'y', label='gaussian, 2')
ax2.plot(x, sf2, 'b', label='splev der=2')
ax2.plot(x, ng2, 'm', label='numpy.gradient 2')
ax2.set_ylabel('second derivative')

plt.title('Original data and 2nd derivatives')
plt.legend()
plt.grid()

无论使用何种方法计算二阶导数，添加到原始函数中的少量噪声都会扰乱二阶导数。