Python中的二阶导数 - scipy/numpy/pandas

Question

我试图在python中使用两个numpy数据数组进行二阶导数.

例如,有问题的数组如下所示:

import numpy as np

x = np.array([ 120. ,  121.5,  122. ,  122.5,  123. ,  123.5,  124. ,  124.5,
        125. ,  125.5,  126. ,  126.5,  127. ,  127.5,  128. ,  128.5,
        129. ,  129.5,  130. ,  130.5,  131. ,  131.5,  132. ,  132.5,
        133. ,  133.5,  134. ,  134.5,  135. ,  135.5,  136. ,  136.5,
        137. ,  137.5,  138. ,  138.5,  139. ,  139.5,  140. ,  140.5,
        141. ,  141.5,  142. ,  142.5,  143. ,  143.5,  144. ,  144.5,
        145. ,  145.5,  146. ,  146.5,  147. ])

y = np.array([  1.25750000e+01,   1.10750000e+01,   1.05750000e+01,
         1.00750000e+01,   9.57500000e+00,   9.07500000e+00,
         8.57500000e+00,   8.07500000e+00,   7.57500000e+00,
         7.07500000e+00,   6.57500000e+00,   6.07500000e+00,
         5.57500000e+00,   5.07500000e+00,   4.57500000e+00,
         4.07500000e+00,   3.57500000e+00,   3.07500000e+00,
         2.60500000e+00,   2.14500000e+00,   1.71000000e+00,
         1.30500000e+00,   9.55000000e-01,   6.65000000e-01,
         4.35000000e-01,   2.70000000e-01,   1.55000000e-01,
         9.00000000e-02,   5.00000000e-02,   2.50000000e-02,
         1.50000000e-02,   1.00000000e-02,   1.00000000e-02,
         1.00000000e-02,   1.00000000e-02,   1.00000000e-02,
         1.00000000e-02,   1.00000000e-02,   5.00000000e-03,
         5.00000000e-03,   5.00000000e-03,   5.00000000e-03,
         5.00000000e-03,   5.00000000e-03,   5.00000000e-03,
         5.00000000e-03,   5.00000000e-03,   5.00000000e-03,
         5.00000000e-03,   5.00000000e-03,   5.00000000e-03,
         5.00000000e-03,   5.00000000e-03])

我现在有f(x) = y,我想要d^2 y / dx^2.

在数值上,我知道我可以插入函数并分析得到导数或使用更高阶的有限差分.我认为有足够的数据可以使用,如果一个或另一个被认为更快,更准确等.

我已经看过np.interp()并且scipy.interpolate没有成功,因为这会返回一个拟合(线性或立方)样条曲线,但不知道如何在此处得到导数.

任何指导都非常感谢.

Answer 1

您可以使用SciPy的的1-d插值数据样条函数.计算样条曲线具有derivative计算导数的便捷方法.

对于您的示例的数据,使用UnivariateSpline给出以下拟合

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

y_spl = UnivariateSpline(x,y,s=0,k=4)

plt.semilogy(x,y,'ro',label = 'data')
x_range = np.linspace(x[0],x[-1],1000)
plt.semilogy(x_range,y_spl(x_range))

适合似乎相当不错,至少在视觉上.您可能想要尝试使用的参数UnivariateSpline.

样条拟合的第二个推导可以简单地获得

y_spl_2d = y_spl.derivative(n=2)

plt.plot(x_range,y_spl_2d(x_range))

结果似乎有点不自然(如果您的数据对应于某些物理过程).您可能想要更改样条拟合参数,改善数据(例如,提供更多样本,执行更少的噪声测量),或者决定分析函数来建模数据并执行曲线拟合(例如,使用sicpy curve_fit)

Answer 2

通过有限差分,y对于数组的每个均值x的一阶导数由下式给出:

dy=np.diff(y,1)
dx=np.diff(x,1)
yfirst=dy/dx

而x的对应值是:

xfirst=0.5*(x[:-1]+x[1:])

对于第二个订单,再次执行相同的过程:

dyfirst=np.diff(yfirst,1)
dxfirst=np.diff(xfirst,1)
ysecond=dyfirst/dxfirst

xsecond=0.5*(xfirst[:-1]+xfirst[1:])