Instantly share code, notes, and snippets.

# pv/splder.py Created May 2, 2013

What would you like to do?
Derivative of B-spline in Python
 r""" Show how to compute a derivative spline. Scipy's splines are represented in terms of the standard B-spline basis functions. In short, a spline of degree k is represented in terms of the knots t and coefficients c by: .. math:: s(x) = \sum_{j=-\infty}^\infty c_{j} B^k_{j}(x) \\ B^0_i(x) = 1, \text{if $t_i \le x < t_{i+1}$, otherwise $0$,} \\ B^k_i(x) = \frac{x - t_i}{t_{i+k} - t_i} B^{k-1}_i(x) + \frac{t_{i+k+1} - x}{t_{i+k+1} - t_{i+1}} B^{k-1}_{i+1}(x) where :math:c_j is nonzero only for 0 <= j <= N and the first and last k knots are at the same position and only used to set up boundary conditions; terms with vanishing denominators are omitted. One can follow standard spline textbooks here, and work out that by differentiating this, one obtains another B-spline, of one order lower: .. math:: s'(x) = \sum_{j=-\infty}^\infty d_j B^{k-1}_j(x) \\ d_j = k \frac{c_{j+1} - c_{j}}{t_{j+k+1} - t_{j+1}} Care needs to be paid at the endpoints: the first knot should be thrown away since the order is reduced by one. """ import numpy as np from scipy.interpolate import UnivariateSpline import matplotlib.pyplot as plt class MyUnivariateSpline(UnivariateSpline): @classmethod def _from_tck(cls, t, c, k): self = cls.__new__(cls) self._eval_args = t, c, k #_data == x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier self._data = [None,None,None,None,None,k,None,len(t),t, c,None,None,None,None] return self def derivative_spline(self): """ Compute the derivative spline of a spline in FITPACK tck representation. """ t, c, k = self._eval_args if k <= 0: raise ValueError("Cannot differentiate order 0 spline") # Compute the denominator in the differentiation formula. dt = t[k+1:-1] - t[1:-k-1] # Compute the new coefficients d = (c[1:-1-k] - c[:-2-k]) * k / dt # Adjust knots t2 = t[1:-1] # Pad coefficient array to same size as knots (FITPACK convention) d = np.r_[d, [0]*k] # Done, return a new spline return self._from_tck(t2, d, k-1) def main(): x = np.linspace(0, 20, 200) y = np.sin(x) # Fit a spline (must use order 4, since roots() works only for order 3) s1 = MyUnivariateSpline(x, y, s=0, k=3+1) # Compute the derivative spline of the fitted spline s2 = s1.derivative_spline() # Roots (works only with order 3 splines) r = s2.roots() print "minima/maxima: x =", r/np.pi print "min/max values: y =", np.sin(r) # Plot and compare the derivative xx = np.linspace(0, 20, 5000) plt.plot(xx, np.cos(xx), 'c-', # exact result xx, s2(xx), 'k--', # derivative spline xx, np.sin(xx), 'm:', # original curve r, s1(r), 'k.', # minima/maxima ) plt.ylim(-1.1, 1.1) plt.savefig('out.png', dpi=96) plt.show() if __name__ == "__main__": main()