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的一维样条函数对数据进行插值。计算的样条曲线具有一种方便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）