[SciPy-User] help interpreting univariate spline
Jason Grout
jason-sage@creativetrax....
Fri Apr 30 07:05:26 CDT 2010
>
> Elliot Hallmark wrote:
>
>> Hi there,
>>
>> I am wanting to use a scipy interpolation routine to generate
>> polynomials approximating function I have sampled. these polynomials
>> will be solved in c (cython), so I need to extract the coefficencts
>> from the interpolation method and pass them to a c function for root
>> finding.
>>
>> using UnivariateSpline I ran this code:
>>
>> import numpy as np
>> from scipy.interpolate import splrep, UnivariateSpline
>>
>> x = np.linspace(1, 10, 100)
>> y = 1/x #approximate a hyperbola
>> g = UnivariateSpline(x, y, w=None, k=3, s=.0005)
>> print g.get_knots()
>> print g.get_coeffs()
>>
>> import matplotlib.pyplot as plt
>> plt.clf()
>> ax = plt.subplot(121)
>> l1, = ax.plot(x, y, "rs", ms=10)
>> l2, = ax.plot(x, g(x), "bo")
>> plt.show()
>>
>>
>>
>>
>> and the output was:
>> [ 1. 2.18181818 3.27272727 5.54545455 10. ]
>>
>> [ 0.98462432 0.66575572 0.423 0.25461626 0.13929847 0.11763985
>> 0.09929961]
>>
>> That is 5 knots and 7 coefficients for a degree=3 spline. naively, I
>> was expecting 4 coefficents for each interval between knots.
>>
>> how can I construct a piecewise polynomial from this output? or can I?
>>
>>
> You might do well to look through the original documentation for Paul
> Dierckx's FITPACK, which scipy.interpolate wraps.
>
>
Also, if you looked at using PPPACK via f2py to get the coefficients of
a cubic spline, you might find the following useful:
http://www.sagenb.org/home/pub/1708/ [1]
Note that in this function, the 4 coefficients returned for each piece
are not the coefficients of x^3, x^2, etc., but rather, the coefficients
are numbers for a Horner-like form of the polynomial (f(x) =
c(1,i)+h*(c(2,i)+h*(c(3,i)+h*c(4,i)/3.)/2.)) where h = x - tau(i) -- see
the documentation for cubspl or the example above to construct the
actual polynomial.)
Thanks,
Jason
[1] (In Sage, using %fortran at the top of a cell wraps the code in an
f2py wrapper.)
More information about the SciPy-User
mailing list