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

Question

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

24

我是一个有用的助手，可以为您翻译文本。

我正在尝试使用两个numpy数据数组在Python中进行二阶导数计算。

例如，所讨论的数组如下：

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，但并没有成功，因为这只返回拟合的（线性或三次）样条曲线，但不知道如何在该点处获得导数。非常感谢您给出的任何指导。

- Jared

1

你看过 np.diff 吗？ - mkhanoyan

我的担忧是我的数据点不是均匀分布的。 - Jared

2个回答

14

通过有限差分，对于您数组中每个x的平均值，y的一阶导数为：

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:])

- manu190466

1

np.diff(np.diff([x*x for x in range(0,10)])) = [2,2,2..]np.diff(np.diff([x*x for x in range(0,10)])) = [2,2,2..] - Blaze

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Stelios · Accepted Answer

你可以使用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）。