Calculating the derivative of points in python

Question

I want to calculate the derivative of points, a few internet posts suggested using np.diff function. However, I tried using np.diff against manually calculated results (chose a random polynomial equation and differentiated it) to see if I end up with the same results. I used the following eq : Y = (X^3) + (X^2) + 7 and the results i ended up with were different. Any ideas why?. Is there any other method to calculate the differential.

In the problem, I am trying to solve, I have recieved data points of fitted spline function ( not the original data that need to be fitted by spline but the points of the already fitted spline). The x-values are at equal intervals. I only have the points and no equation, what I require is to calculate, the first, second and third derivatives. i.e dy/dx, d2y/dx2, d3y/dx3. Any ideas on how to do this?. Thanks in advance.

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]

The values of numpy diff method [10,24,44,60] is different from yval_dashList = [5,16,35,64,107]

Answer 1

I like @lgsp's answer. I will add that you can directly estimate the derivative without having to worry about how much space there is between the values. This just uses the symmetric formula for calculating finite differences, described at this wikipedia page.

Take note, though, of the way delta is specified. I found that when it is too small, higher-order estimates fail. There's probably not a 100% generic value that will always work well!

Also, I simplified your code by taking advantage of numpy broadcasting over arrays to eliminate for loops.

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

And the outputs:

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]

Answer 2

The idea behind what you are trying to do is correct, but there are a couple of points to make it work as intended:

1. There is a typo in calc_diffY(X), the derivative of X**2 is 2*X, not 2**X:

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

By doing this you don't obtain much better results:

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

2. To calculate the numerical derivative you should do a "Difference quotient" which is an approximation of a derivative

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

The approximation becomes better and better if the values of the points are more dense. The difference between your points on the x axis is 1, so you end up in this situation (in blue the analytical derivative, in red the numerical):

If you reduce the difference in your x points to 0.1, you get this, which is much better:

Just to add something to this, have a look at this image showing the effect of reducing the distance of the points on which the derivative is numerically calculated, taken from Wikipedia: