[SciPy-User] help interpreting univariate spline

Jason Grout jason-sage at creativetrax.com
Fri Apr 30 08:05:26 EDT 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