[SciPy-User] Newbie Question :: Natural Cubic Spline

Etrade Griffiths etrade.griffiths@dsl.pipex....
Thu Oct 22 03:17:16 CDT 2009


>Hello:
>
>I am new to Python and the SciPy libraries and I was wondering if 
>someone could help me with the following.
>
>Given the following x and y : values
>
>x_list = 
>[33,              56,     56.00000000002,       147,    238,    329, 
>    420,    511,    602,    693,    791]
>
>y_list = 
>[0.99974,    0.99949,            0.99949, 
>0.99816,0.99631,0.99383,0.99043,0.98610,0.98078,0.97460,0.96704]
>
>
>I want to use a Natural Cubic Spline in order to determine the 
>points at the following values:
>
>interp_x_points    = 
>[31,62,92,123,153,184,215,243,274,304,335,365,396,427,457,488,518,549,580,608,639,669,700] 
>
>
>I have looked through the documentation and have had a tough time 
>determining the correct syntax or which function to use.
>
>I was wondering what would be the SciPy syntax for the following:
>Also, does scipy have the ability to extrapolate if a give X value 
>is outside a specified range?
>In this example, please notice interp_x_poiints has a value (31) 
>which is outside of the x_list range.

Try this:

# Example use of cubic spline

import scipy.interpolate

x_list = [33, 56, 56.00000000002, 147, 238, 329, 420, 511, 602, 693, 791]
y_list = [0.99974, 0.99949, 
0.99949,0.99816,0.99631,0.99383,0.99043,0.98610,0.98078,0.97460,0.96704]

interp_x_points = 
[31,62,92,123,153,184,215,243,274,304,335,365,396,427,457,488,518,549,580,608,639,669,700]

# Get the knot points (s=0.0 => "natural" splines)

tck = scipy.interpolate.splrep(x_list, y_list, s=0.0)

for xval in interp_x_points:
     try:
         yval = scipy.interpolate.splev(xval,tck,der=0)
     except:
         yval = -1.0E30

     print "x value: %5.1f; y value:  %7.5f" % (xval, yval)

Spline fitting is a two step process: first calculate the knot points 
and then do the interpolation.  Press et al have a nice description 
of the process here (Chapter 3.3)

http://www.nrbook.com/a/bookcpdf.php

Note that you may have to install the free plug in mentioned on this web page

HTH

Alun Griffiths










More information about the SciPy-User mailing list