Derivatives blow up in python

Question

I am trying to find higher order derivatives of a dataset (x,y). x and y are 1D arrays of length N.

Let's say I generate them as :

xder0=np.linspace(0,10,1000)
yder0=np.sin(xder0)

I define the derivative function which takes in 2 array (x,y) and returns (x1, y1) where y1 is the derivative calculated at each index as : (y[i+1]-y[i])/(x[i+1]-x[i]). x1 is just the mean of x[i+1] and x[i]

Here is the function that does it:

def deriv(x,y):
    delx =np.zeros((len(x)-1), dtype=np.longdouble)
    ydiff=np.zeros((len(x)-1), dtype=np.longdouble)
    for i in range(len(x)-1):
            delx[i]  =(x[i+1]+x[i])/2.0
            ydiff[i] =(y[i+1]-y[i])/(x[i+1]-x[i])
    return delx, ydiff

Now to calculate the first derivative, I call this function as:

xder1, yder1 = deriv(xder0, yder0)

Similarly for second derivative, I call this function giving first derivatives as input:

xder2, yder2 = deriv(xder1, yder1)

And it goes on:

xder3, yder3 = deriv(xder2, yder2)
xder4, yder4 = deriv(xder3, yder3)
xder5, yder5 = deriv(xder4, yder4)
xder6, yder6 = deriv(xder5, yder5)
xder7, yder7 = deriv(xder6, yder6)
xder8, yder8 = deriv(xder7, yder7)
xder9, yder9 = deriv(xder8, yder8)

Something peculiar happens after I reach order 7. The 7th order becomes very noisy! Earlier derivatives are all either sine or cos functions as expected. However 7th order is a noisy sine. And hence all derivatives after that blow up.

Any idea what is going on?

Answer 1

This is a well known stability issue with numerical interpolation using equally-spaced points. Read the answers at http://math.stackexchange.com.

To overcome this problem you have to use non-equally-spaced points, like the roots of Lagendre polynomial. The instability occurs due to the unavailability of information at the boundaries, thus more concentration of points at the boundaries is required, as per the roots of say Lagendre polynomials or others with similar properties such as Chebyshev polynomial.