[SciPy-user] cubic splines in scipy
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
>(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