Мне нравится ответ @ lgsp. Я добавлю, что вы можете напрямую оценить производную, не беспокоясь о том, сколько места между значениями. Это просто использует формулу симметрии c для расчета конечных разностей, описанную в на этой странице википедии .
Обратите внимание, что указан способ delta
. Я обнаружил, что когда он слишком мал, оценки более высокого порядка терпят неудачу. Вероятно, нет 100% общего значения c, которое всегда будет работать хорошо!
Кроме того, я упростил ваш код, воспользовавшись преимуществом 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]