# public pv / splder.py Created 2013-05-02

Derivative of B-spline in Python

splder.py
Python
 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 r"""Show how to compute a derivative spline. Scipy's splines are represented in terms of the standard B-spline basisfunctions. In short, a spline of degree k is represented in terms of theknots 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 firstand last k knots are at the same position and only used to set upboundary conditions; terms with vanishing denominators are omitted. One can follow standard spline textbooks here, and work out that bydifferentiating this, one obtains another B-spline, of one orderlower: .. 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 knotshould be thrown away since the order is reduced by one. """import numpy as npfrom 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()