在python中计算点的导数

Question

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

在我试图解决的问题中，我收到了拟合样条函数的数据点（不是需要用样条拟合的原始数据，而是已经拟合的样条函数的点）。x值的间隔相等。我只有点，没有方程，我需要的是计算，一阶，二阶和三阶导数。i、 e 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 method[10,24,44,60]的值不同于yval_dashList = [5,16,35,64,107]

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.]

1. 要计算数值导数，你应该做一个“差商”，它是导数的近似值

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

如果点的值更密集，则近似效果会越来越好。 x轴上的点之间的差值为1，因此在这种情况下结束（蓝色表示分析导数，红色表示数值）：

如果您将x点的差值减少到0.1，您会得到这样的结果，这要好得多：

为了增加一些内容，请看这张图片，它显示了减少数值计算导数的点的距离的效果，取自Wikipedia：

Answer 2

我喜欢@lgsp的答案。我要补充的是，你可以直接估计导数，而不必担心值之间有多少空间。这只是使用对称公式来计算有限差分，如this wikipedia page所述

不过，请注意delta的指定方式。我发现当它太小时，高阶估计会失败。可能没有一个100%的通用值能够始终很好地工作

此外，我还利用数组上的numpy广播来消除for循环，从而简化了代码

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]