在Python中计算点的导数

Question

我想计算点的导数，一些互联网文章建议使用np.diff函数。然而，我尝试使用np.diff对手动计算的结果进行比较（选择一个随机的多项式方程并对其进行微分），以查看是否得到相同的结果。我使用以下方程：Y =（X ^ 3）+（X ^ 2）+ 7，但我最终得到的结果不同。有什么想法吗？是否有其他方法来计算微分值。

在我要解决的问题中，我收到了拟合样条函数的数据点（不是需要通过样条拟合的原始数据，而是已经拟合的样条的点）。x值处于等间隔。我只有这些点，没有方程式，我需要计算一阶、二阶和三阶导数。即dy/dx，d2y/dx2，d3y/dx3。有关如何执行此操作的任何想法吗？先谢谢您。

xval = [1,2,3,4,5]
yval = []
yval_dashList = []

#selected a polynomial equation
def calc_Y(X):
      Y = (X**3) + (X**2) + 7
      return(Y)

#calculate y values using equatuion 
for i in xval:
    yval.append(calc_Y(i))

#output: yval = [9,19,43,87,157]

#manually differentiated the equation or use sympy library (sym.diff(x**3 + x**2 + 7))
def calc_diffY(X):
   yval_dash = 3*(X**2) + 2**X

#store differentiated y-values in a list
for i in xval:
    yval_dashList.append(yval_dash(i))

#output: yval_dashList = [5,16,35,64,107]

#use numpy diff method on the y values(yval)
numpyDiff = np.diff(yval)
#output: [10,24,44,60]

numpy的diff方法计算[10,24,44,60]的差异，与yval_dashList = [5,16,35,64,107]不同。

- imantha

---

可能的答案是 numpy.diff 不会做你期望的事情。numpy.diff 仅仅告诉你数组中相邻数值的差异。在一些情况下，这可能近似于函数的导数，但大多数情况下并不是这样。 - senderle

那么...你的目标实际上是基于函数本身来数值估计导数吗？ - senderle

你可以使用 autograd 模块。 - hilberts_drinking_problem

Answer 1

你尝试做的想法是正确的，但有几点需要注意才能实现预期效果：

1. calc_diffY（X）中有一个拼写错误，X ** 2的导数是2 * X，而不是2 ** X：

  def calc_diffY(X):    
      yval_dash = 3*(X**2) + 2*X

通过这样做，你不会得到更好的结果：

yval_dash = [5, 16, 33, 56, 85]
numpyDiff = [10. 24. 44. 70.]

要计算数值导数，您应该使用“差商”，它是导数的近似值。

numpyDiff = np.diff(yval)/np.diff(xval)

逼近会变得越来越好，如果点的值更加密集。您在x轴上的点之间的差为1，因此您会遇到这种情况（蓝色是解析导数，红色是数值导数）：

如果你将x点之间的差异减少到0.1，你会得到更好的结果：

补充一点，看看这张图片，它展示了从Wikipedia中获取的通过减小数值求导计算点之间距离的影响：

- lgsp

Answer 2

我喜欢@lgsp的答案。此外，您可以直接估计导数而无需担心值之间有多少间隔。这只是使用对称公式来计算有限差分，在这个维基百科页面中有描述。

请注意，delta的指定方式。我发现当它太小时，高阶估计就会失败。可能没有一个100％通用的值总是表现良好！

此外，我通过利用numpy广播数组来消除循环，简化了您的代码。

import numpy as np

# selecte a polynomial equation
def f(x):
    y = x**3 + x**2 + 7
    return y

# manually differentiate the equation
def f_prime(x):
    return 3*x**2 + 2*x

# numerically estimate the first three derivatives
def d1(f, x, delta=1e-10):
    return (f(x + delta) - f(x - delta)) / (2 * delta)

def d2(f, x, delta=1e-5):
    return (d1(f, x + delta, delta) - d1(f, x - delta, delta)) / (2 * delta)

def d3(f, x, delta=1e-2):
    return (d2(f, x + delta, delta) - d2(f, x - delta, delta)) / (2 * delta)

# demo output
# note that functions operate in parallel on numpy arrays -- no for loops!
xval = np.array([1,2,3,4,5])

print('y  = ', f(xval))
print('y\' = ', f_prime(xval))
print('d1 = ', d1(f, xval))
print('d2 = ', d2(f, xval))
print('d3 = ', d3(f, xval))

输出结果如下：

y  =  [  9  19  43  87 157]
y' =  [ 5 16 33 56 85]
d1 =  [ 5.00000041 16.00000132 33.00002049 56.00000463 84.99995374]
d2 =  [ 8.0000051  14.00000116 20.00000165 25.99996662 32.00000265]
d3 =  [6.         6.         6.         6.         5.99999999]

- senderle

---

一种看待这个问题的方式是注意到，如果对结果考虑的x位置是计算差异的平均点，则由@AccLok计算的导数看起来很好，也就是说：xval=[1.5, 2.5, 3.5, 4.5]numpyDiff = [10, 24, 44, 70,] - lgsp

@lgsp 抱歉，刚回来看到这个消息。这是一个很好的观点 - 实际上，有一个简单的变化可以很好地解决这个问题：(np.diff(yval[1:]) + np.diff(yval[:-1])) / 2)。它只给出中间三个值的结果，但每次只差1。 - senderle