Конечные разностные функции на языке Python?

Я просматривал Numpy/Scipy для модулей, содержащих конечные разностные функции. Тем не менее, ближайшая вещь, которую я обнаружил, numpy.gradient(), которая хороша для конечных разностей 1-го порядка точности 2-го порядка, но не так много, если вам нужны производные более высокого порядка или более точные методы. Я даже не нашел очень много специальных модулей для такого рода вещей; большинство людей, похоже, делают "рулонные", поскольку они в них нуждаются. Поэтому мой вопрос заключается в том, что кто-нибудь знает какие-либо модули (либо часть Numpy/Scipy, либо сторонний модуль), которые специально посвящены методам конечных разностей с более высоким порядком (как по точности, так и по производным). У меня есть собственный код, над которым я работаю, но в настоящее время он медленный, и я не буду пытаться его оптимизировать, если там что-то уже доступно.

Обратите внимание, что я говорю о конечных разностях, а не о производных. Я видел как scipy.misc.derivative(), так и Numdifftools, которые берут производную от аналитической функции, которой у меня нет.

Ответ 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)

Ответ 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

Ответ 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()

Обратите внимание на нулевые пересечения первой производной (красная кривая) на вершинах и впадинах сплайна (черная кривая).