2
$\begingroup$

Question

Below is a plot of a graph $y$ and its derivative $dy/dx$ calculated using python's numpy.gradientwhich approximates the derivative with finite differences using a central scheme. Clearly there is a problem in that unphysical kinks have been introduced in the result.

Since differentiation will act to amplify any high frequency content in the original data, numerical differentiation is always tricky although here my data is relatively smooth and noise free. Does anyone recognise this type of behaviour and what can be done to get a proper derivative?

Some info on the data:

$y$ is data that is sampled from a computational mesh and has many more sample points for small $x$, getting coarser as $x$ increases. The data you see plotted is found by linearly interpolating on a uniform number of data points. However because of the non-uniform spacing I have had to use a large number of points so that the data points for small $x$ are effectively the "true" data values, but at high $x$ there are many points in between the "true" ones that are just linear interpolations, adding no real information to the graph. At first I thought maybe there is some floating point rounding problems that are causing

enter image description here

EDIT: MCVE

To illustrate the problem here is a small subset of the data to produce this figure:

enter image description here

import numpy as np
import matplotlib.pyplot as plt

# Subset of data such that
# sub_data = data[50:150,:]

sub_data = np.asarray([[1.700089160294737667e-05, 9.632603526115417480e-01],
[1.733388671709690243e-05, 9.813255667686462402e-01],
[1.766688365023583174e-05, 9.993913173675537109e-01],
[1.799987876438535750e-05, 1.017457604408264160e+00],
[1.833267015172168612e-05, 1.035513401031494141e+00],
[1.866566708486061543e-05, 1.053580760955810547e+00],
[1.899866219901014119e-05, 1.071648478507995605e+00],
[1.933165913214907050e-05, 1.089716911315917969e+00],
[1.966444870049599558e-05, 1.107774615287780762e+00],
[2.000245694944169372e-05, 1.126114964485168457e+00],
[2.033545388258062303e-05, 1.144184470176696777e+00],
[2.066844899673014879e-05, 1.162254691123962402e+00],
[2.100123856507707387e-05, 1.180314183235168457e+00],
[2.133423549821600318e-05, 1.198385119438171387e+00],
[2.166586818930227309e-05, 1.216382265090942383e+00],
[2.199886512244120240e-05, 1.234277606010437012e+00],
[2.233165650977753103e-05, 1.251089692115783691e+00],
[2.266465162392705679e-05, 1.267913579940795898e+00],
[2.299764855706598610e-05, 1.284739017486572266e+00],
[2.333043812541291118e-05, 1.301555633544921875e+00],
[2.366844637435860932e-05, 1.318636536598205566e+00],
[2.400144148850813508e-05, 1.335466384887695312e+00],
[2.433443842164706439e-05, 1.352297544479370117e+00],
[2.466722798999398947e-05, 1.369119763374328613e+00],
[2.500022492313291878e-05, 1.385953545570373535e+00],
[2.533322003728244454e-05, 1.402788400650024414e+00],
[2.566621697042137384e-05, 1.419624686241149902e+00],
[2.599900835775770247e-05, 1.436451792716979980e+00],
[2.633200347190722823e-05, 1.453290104866027832e+00],
[2.666500040504615754e-05, 1.470129728317260742e+00],
[2.699799551919568330e-05, 1.486970305442810059e+00],
[2.733579822233878076e-05, 1.504054188728332520e+00],
[2.766743091342505068e-05, 1.520827651023864746e+00],
[2.800042602757457644e-05, 1.537671327590942383e+00],
[2.833321741491090506e-05, 1.554505705833435059e+00],
[2.866621252906043082e-05, 1.571351289749145508e+00],
[2.899920946219936013e-05, 1.588197588920593262e+00],
[2.933220639533828944e-05, 1.605044960975646973e+00],
[2.966499596368521452e-05, 1.621882915496826172e+00],
[2.999799289682414383e-05, 1.638731956481933594e+00],
[3.033098801097366959e-05, 1.655581831932067871e+00],
[3.066398494411259890e-05, 1.672432422637939453e+00],
[3.100178582826629281e-05, 1.689526319503784180e+00],
[3.133478094241581857e-05, 1.706378579139709473e+00],
[3.166777969454415143e-05, 1.723231554031372070e+00],
[3.200077480869367719e-05, 1.740085244178771973e+00],
[3.233356619603000581e-05, 1.756929397583007812e+00],
[3.266656131017953157e-05, 1.773784637451171875e+00],
[3.299955642432905734e-05, 1.790640354156494141e+00],
[3.333234781166538596e-05, 1.806897401809692383e+00],
[3.366534292581491172e-05, 1.822490572929382324e+00],
[3.399697743589058518e-05, 1.838021636009216309e+00],
[3.433498568483628333e-05, 1.853852272033691406e+00],
[3.466777343419380486e-05, 1.869441509246826172e+00],
[3.500077218632213771e-05, 1.885041832923889160e+00],
[3.533376730047166348e-05, 1.900644183158874512e+00],
[3.566676605259999633e-05, 1.916247844696044922e+00],
[3.599955380195751786e-05, 1.931843996047973633e+00],
[3.633255255408585072e-05, 1.947450876235961914e+00],
[3.666554766823537648e-05, 1.963059425354003906e+00],
[3.699854278238490224e-05, 1.978669643402099609e+00],
[3.733133416972123086e-05, 1.994271874427795410e+00],
[3.766432928387075663e-05, 2.009885072708129883e+00],
[3.800233753281645477e-05, 2.025733470916748047e+00],
[3.833533264696598053e-05, 2.041349411010742188e+00],
[3.866812403430230916e-05, 2.056957483291625977e+00],
[3.900111914845183492e-05, 2.072576045989990234e+00],
[3.933411426260136068e-05, 2.088196039199829102e+00],
[3.966690564993768930e-05, 2.103807926177978516e+00],
[3.999854016001336277e-05, 2.119366645812988281e+00],
[4.033153527416288853e-05, 2.134990215301513672e+00],
[4.066453038831241429e-05, 2.150615453720092773e+00],
[4.099732177564874291e-05, 2.166232347488403320e+00],
[4.133032052777707577e-05, 2.181859731674194336e+00],
[4.166832513874396682e-05, 2.197722673416137695e+00],
[4.200132389087229967e-05, 2.213352441787719727e+00],
[4.233411164022982121e-05, 2.228973865509033203e+00],
[4.266711039235815406e-05, 2.244606018066406250e+00],
[4.300010550650767982e-05, 2.260239124298095703e+00],
[4.333310062065720558e-05, 2.275873422622680664e+00],
[4.366589200799353421e-05, 2.291499376296997070e+00],
[4.399888712214305997e-05, 2.307135581970214844e+00],
[4.433188223629258573e-05, 2.322772979736328125e+00],
[4.466467362362891436e-05, 2.338401794433593750e+00],
[4.499766873777844012e-05, 2.352876663208007812e+00],
[4.533567698672413826e-05, 2.367459297180175781e+00],
[4.566867210087366402e-05, 2.381829023361206055e+00],
[4.600146712618879974e-05, 2.396192073822021484e+00],
[4.633309799828566611e-05, 2.410506963729858398e+00],
[4.666609311243519187e-05, 2.424882888793945312e+00],
[4.699908822658471763e-05, 2.439260959625244141e+00],
[4.733187961392104626e-05, 2.453632116317749023e+00],
[4.766487472807057202e-05, 2.468014001846313477e+00],
[4.799787348019890487e-05, 2.482398033142089844e+00],
[4.833086859434843063e-05, 2.496783494949340820e+00],
[4.866365998168475926e-05, 2.511162519454956055e+00],
[4.900166459265165031e-05, 2.525767564773559570e+00],
[4.933466334477998316e-05, 2.540158748626708984e+00],
[4.966766209690831602e-05, 2.554551839828491211e+00],
[5.000044984626583755e-05, 2.568938016891479492e+00]])


x = sub_data[:,0]
y = sub_data[:,1]
dx = x[1]-x[0]

dydx = np.gradient(y,dx)

plt.figure()
fig, ax1 = plt.subplots(figsize=(8,4))
ax2 = ax1.twinx()

ax1.plot(x,y,'k-')
ax2.plot(x,dydx, 'r-')
plt.show()
$\endgroup$
  • $\begingroup$ The jaggedness in that first derivative, when looking at the plot for $Y$, isn't something I would expect. Have you tried computing a numerical derivative yourself and seeing what you get? $\endgroup$ – spektr Sep 28 '16 at 16:43
  • $\begingroup$ @choward I have tried using the simplest thing I could think of which is forward differencing (calculating the difference array of y and dividing it by the difference array of x) but the result is still the same! Im thinking the problem is in my data and the fact that I have so many closely sampled points...? $\endgroup$ – Dipole Sep 28 '16 at 16:46
  • $\begingroup$ @Jack You could fit your data $y(x)$ with a differentiable function (e.g., a polynomial) in order to find a better estimate of the derivative. Scipy's scipy.interpolate.UnivariateSpline provides such a fit and has a convenient built in derivative() method $\endgroup$ – Stelios Sep 28 '16 at 16:50
  • 1
    $\begingroup$ Can you post an MCVE? $\endgroup$ – Kirill Sep 28 '16 at 16:51
  • $\begingroup$ @Stelios I actually tried that already, but that spline returns all NaNs (I have also removed any possible duplicate points already) $\endgroup$ – Dipole Sep 28 '16 at 16:52
6
$\begingroup$

The $x$-coordinates of your data points are not equally spaced (x[1:]-x[:-1] is not constant), so numpy.gradient is not applicable because it assumes that the data is equally spaced (http://docs.scipy.org/doc/numpy/reference/generated/numpy.gradient.html). Even just forward differences on their own would be better than using an inaccurate value of $\Delta x$ in centered differences.

With forward differences, the small bumps disappear, but the discontinuities remain, so those probably come from the data itself. You can smooth them by constructing a spline interpolant (http://docs.scipy.org/doc/scipy/reference/interpolate.html).

$\endgroup$
  • $\begingroup$ Darn! Sadly I knew this but I had assumed that due to the way I had acquired the data the points were equally spaced...Thanks! $\endgroup$ – Dipole Sep 28 '16 at 20:12
2
$\begingroup$

Here's a plot of the derivative evaluated by first performing an spline fit of the data. Compared to the finite difference approach, the result appears more natural.

from scipy.interpolate import UnivariateSpline

spl = UnivariateSpline(sub_data[:,0],sub_data[:,1])

x_range = np.linspace(sub_data[0,0], sub_data[-1,0],1000)
plt.plot(x_range,spl.derivative(1)(x_range),'b', label = 'spline')
plt.plot(x,dydx, 'r-', label = 'finite diff.')
plt.legend(loc = 'best')

enter image description here

$\endgroup$
  • $\begingroup$ Thanks - I had originally used splev without success for some reason but this worked nicely on my full data although I had to play around with the smoothing factor a bit! $\endgroup$ – Dipole Sep 28 '16 at 20:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.