[SciPy-user] cubic splines in scipy

Travis Oliphant oliphant.travis at ieee.org
Sat Jan 28 18:06:21 CST 2006

Andrew Straw wrote:

>(Warning: non math guru alert.)
>I'm looking to do 1D interpolation using cubic splines:
>http://mathworld.wolfram.com/CubicSpline.html That page has a 2D
>picture, but the math is for 1D. The important bit for me is that
>interpolated curve passes through the control points and the whole thing
>is continuously differentiable.
There are two incarnations of splines in SciPy.  The first is from 
fitpack.  It is the traditional way to view cubic splines as it allows 
for control-points to be placed in arbitrary ways.

The specialized cubic-spline algorithms in scipy.signal make use of 
speed enhancements that are possible when your control points are 
*equally-spaced* and you are willing to assume specific kinds of 
boundary conditions (Michael Unser has written extensively about this 
use of cubic splines and his papers form the basis for what is in 
scipy.signal).   The big win is that by making these assumptions your 
matrix-to-be inverted is circulant and the inversion can be accomplished 
very quickly.

>(For what it's worth, I need to interpolate discretely sampled data from
>a discrete-time simulation that I then want to use as forcing functions
>for use in a scipy ODE solver. So, my discretely sampled data contain no
>noise and I don't want to smooth them, I'm just looking for intermediate
>values that are continuously differentiable.)
Look in the old scipy tutorial for discussion on the spline functions.  
Both, the interpolate in fitpack and the B-spline approach are 
discussed.  If you have equally-spaced points (it looks like you do), 
you can use the B-spline approach.  Do not use the Mathworld discussion 
for this approach (it's related but the formulas there won't necessarily 

What to do instead?  Well, the bsplines need some more helper functions 
to make what you want to do easy. 

Instead, use the fact that the acutal interpolating function is:

y(x) = sum(c_j*beta_3(x/deltax - j))  

where beta_3 is a cubic-spline interpolator function (a symmetric 
kernel)  which scipy.signal.bspline(x/deltax-j,3)  will evaluate for 
you, and cj are the (spline coefficients) returned from cspline1d(y). 

In order to get the edges to come out right you will need to extend the 
returned coefficients using mirror symmetric boundary conditions and 
include those extra knots in your sum near the edges.

We really need a faster way to compute the interpolation function 
(besides evaluating beta_3(x/deltax-j) at every new x and j point and 
summing the results because most of them are zero).    This would be a 
good addition and a quick-way to get 2-d interpolation working quickly.


More information about the SciPy-user mailing list