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数组,那么以下内容适用于:scipy.ndimage如果你想为ndarrays做这个,它就内置在包中.但要小心,因为当然这不会给你完整的渐变但我相信所有方向的产物.有更好专业知识的人希望能说出来.

这是一个例子:

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

绝对像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]')

... 像这样使用'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')

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

Answer 3

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

• 沿任意轴微分任意维数的数组
• 任何所需顺序的偏导数
• 来自矢量微积分的标准算子，如梯度、散度和卷曲
• 可以处理均匀和非均匀网格
• 可以处理具有常数和可变系数的导数的任意线性组合
• 可以指定精度顺序
• 速度完全矢量化
• 计算均匀和非均匀网格的任何阶次和精度的原始有限差分系数

Answer 4

另一种方法是区分数据的插值。这是 unutbu 建议的，但我没有看到此处或任何链接问题中使用的方法。UnivariateSpline例如， fromscipy.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()

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