Search code examples
pythonnumpyscipyderivative

Get derivative of data in python

I write a program to get derivative. InterpolatedUnivariateSpline is used for calculating f(x+h). The red line is derivative of cosine, the green line is cosine consine, the blue line is -sine function. The red and blue line are matched. It works well in the following.

Result

from scipy.interpolate import InterpolatedUnivariateSpline
import numpy as np 
import matplotlib.pyplot as plt

pi = np.pi
x = np.arange(0,5*pi,0.2*pi)
y = np.cos(x) 
f2 = InterpolatedUnivariateSpline(x, y)
#Get dervative
der = []
for i in range(len(y)):

    h = 1e-4
    der.append( ( f2(x[i]+h)-f2(x[i]-h) )/(2*h) )
der = np.array(der)   

plt.plot(x, der, 'r', x, y, 'g', x, -np.sin(x),'b')
plt.show()

But I encounter some problem. In my project, my variable x(frequency) vary from 10^7 to 2.2812375*10^9, its step is 22487500,so I change array x. As a result, I get the following result.

Result2

The derivative is a line and almost close to 0, It is not -sine function. How do I solve this?

Solution

  • You have a loss of significance problem. It means that when adding a large floating point number to a small one, the precision of the small one is partially lost as a numpy double can only hold 64 bits of information.

    To solve this issue you have to make sure that the scale of the numbers you add/multiply/divide is not too different. One simple solution is dividing x by 1e9 or multiplying h by 1e9. If you do this you get essentially the same precision as in your example.

    Also if x has high enough resolution a simpler way to numerically differentiate the function would be der = np.diff(y) / np.diff(x). This way you do not have to worry about setting h. However in this case note that dy is one element shorter that y, and dy[i] is actually an approximation of the derivative at `(x[i] + x[i+1]) / 2. So to plot it you would do:

    der = np.diff(y) / np.diff(x)
x2 = (x[:-1] + x[1:]) / 2
plt.plot(x2, der, 'r', x, y, 'g', x, -np.sin(x),'b')