<p>我喜欢@lgsp的答案。我要补充的是,你可以直接估计导数,而不必担心值之间有多少空间。这只是使用对称公式来计算有限差分,如<a href="https://en.wikipedia.org/wiki/Numerical_differentiation" rel="noreferrer">this wikipedia page</a>所述</p>
<p>不过,请注意<code>delta</code>的指定方式。我发现当它太小时,高阶估计会失败。可能没有一个100%的通用值能够始终很好地工作</p>
<p>此外,我还利用数组上的<code>numpy</code>广播来消除for循环,从而简化了代码</p>
<pre><code>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))
</code></pre>
<p>以及产出:</p>
<pre><code>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]
</code></pre>