Python有限差分函数？答案

【问题标题】：Python finite difference functions?Python有限差分函数？
【发布时间】：2013-09-30 05:56:29
【问题描述】：

我一直在 Numpy/Scipy 中寻找包含有限差分函数的模块。但是，我发现的最接近的东西是numpy.gradient()，它适用于二阶精度的一阶有限差分，但如果您想要高阶导数或更精确的方法，则效果不佳。我什至还没有为这类事情找到很多特定的模块。大多数人似乎都在根据需要做“自己动手”的事情。所以我的问题是，是否有人知道专门用于高阶（精度和导数）有限差分方法的任何模块（Numpy/Scipy 的一部分或第三方模块）。我有自己的代码正在处理，但它目前有点慢，如果已经有可用的东西，我不会尝试对其进行优化。

请注意，我说的是有限差分，而不是导数。我见过scipy.misc.derivative() 和Numdifftools，它们采用分析函数的导数，而我没有。

【问题讨论】：

对数据进行插值，然后在插值函数上使用scipy.misc.derivative 对您有用吗？
看看numpy.diff (docs.scipy.org/doc/numpy/reference/generated/numpy.diff.html)。
另外，如果您的数据是周期性的，则有scipy.fftpack.diff (docs.scipy.org/doc/scipy/reference/generated/…)。这里有一个使用示例：wiki.scipy.org/Cookbook/KdV
@unutbu：我想理论上它可以工作，如果我正确阅读文档，通过将 dx 参数设置为与网格间距相同并强制 x0 位于其中之一网格点。但是，我的意图是将其用于 2 维和 3 维数组，因此对于 MxNxP 数组，我必须为 M 方向的导数创建 NxP interp1d 对象，这似乎有点慢。它还需要一个 Python 循环，出于同样的原因，我想避免这种循环。
@askewchan：这里的精度是指由于泰勒级数截断而导致的有限差分近似中的误差。更精确的有限差分方法保留更多的泰勒级数项，因此更接近该点的真实导数。一阶保留的术语少于二阶，依此类推。

标签： python numpy scipy

【解决方案1】：

FastFD 可以提供帮助：

pip install fastfd

基本文档可在此处找到：

https://github.com/stefanmeili/FastFD

这里：

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

【讨论】：

【解决方案2】：

快速做到这一点的一种方法是使用高斯核的导数进行卷积。简单的情况是你的数组与[-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)

【讨论】：

【解决方案3】：

您可能想看看findiff project。我自己试过了，它让您可以方便地获取任何维度、任何导数顺序和任何所需精度顺序的 numpy 数组的导数。该项目网站说它具有：

沿任意轴区分任意维数的数组
任意阶的偏导数
矢量演算中的标准运算符，如梯度、散度和旋度
可以处理均匀和不均匀的网格
可以处理具有恒定和可变系数的导数的任意线性组合
可以指定精度顺序
完全矢量化以提高速度
为均匀和非均匀网格计算任意阶和精度的原始有限差分系数

【讨论】：

【解决方案4】：

另一种方法是区分数据的插值。这是 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()

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

【讨论】：

请记住，UnivariateSpline 仅限于五度样条 (k<=5)，这意味着您仅限于 n=4 导数。此外，与高斯卷积一样，需要根据需要更改平滑因子。

【解决方案5】：

绝对喜欢 askewchan 给出的答案。这是一项很棒的技术。但是，如果您需要使用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]')

...在numpy.convolve 中使用'same' 选项，如下所示：

#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-index。 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')

使用 s/bot/by/g 编辑更正的自动完成

【讨论】：