Python有限差分函数？

Question

我在Numpy/Scipy中四处寻找包含有限差分函数的模块。然而，我发现的最接近的东西是numpy.gradient()，它有利于二阶精度的一阶有限差分，但如果你想要更高的导数或更精确的方法，那就不多了。我甚至没有为这类事情找到很多特定的模块；大多数人似乎都在做一件“自己拥有”的事情，因为他们需要它们。因此，我的问题是，是否有人知道任何模块( Numpy/Scipy或第三方模块的一部分)专门用于高阶(无论是精度还是导数)有限差分方法。我有自己的代码，我正在做，但它目前有点慢，我不会试图优化它，如果有什么东西已经可用。

注意，我说的是有限差分，而不是导数。我看到了scipy.misc.derivative()和数扩散工具，它们接受分析函数的导数，但我没有。

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

绝对喜欢阿斯凯文的回答。这是一种很棒的技术。但是，如果您需要使用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]')

...use '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]')

...to避免过一次索引错误.

在绘图时也要小心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')

​

​

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

Answer 3

您可能想看看findiff项目。我自己也试过了，让我们方便地获取任意维数的numpy数组的导数，任何导数顺序，以及任何期望的精度顺序。该项目的网站说，它的特点是：

• 沿任意轴区分任意维数的数组
• 任意阶偏导数
• 向量演算中的标准算子，如梯度、散度和卷曲
• 能处理均匀和非均匀网格。
• 可以处理常系数和变系数导数的任意线性组合。
• 可以指定精度顺序。
• 完全矢量化的速度
• 计算均匀和非均匀网格的任意阶和精度的原始有限差分系数。