Skip to content

Instantly share code, notes, and snippets.

@pv
Created May 2, 2013 18:43
  • Star 4 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save pv/5504366 to your computer and use it in GitHub Desktop.
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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment