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?

The difference formulae for approximating derivations are all *ill-conditioned*, because we are substracting two relatively large (`x[i]`, `x[i+1]`, `y[i]`, `y[i + 1]`) numbers multiple times to get a relatively small difference. `order = 7` seems to hit the dropoff-point for you. — dhke, Dec 03 '15 at 09:17
Look for a library that does handles the numerical part. [numdifftools](https://github.com/pbrod/numdifftools) comes to mind. As you can see, a direct, unstabilized approach does only get you so far. — dhke, Dec 03 '15 at 10:08
One possibility would be to fit some function to your data that you can then take analytical derivatives of. — ali_m, Dec 03 '15 at 13:31
I would recommend reading Section 5.7 of [Numerical Recipes](http://apps.nrbook.com/empanel/index.html#). It gives multiple methods of computing numerical derivatives while minimizing the stability issues. You have access to a couple of pages per month. Should be enough for your purposes, I would think. — Joey Dumont, Dec 04 '15 at 03:00

score 2 · Answer 1 · edited Apr 13 '17 at 12:19

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.

Derivatives blow up in python

1 Answers1