[SciPy-user] SciPy cubic interpolation coefficients
josef.pktd@gmai...
josef.pktd@gmai...
Mon Jun 22 14:07:30 CDT 2009
On Sun, Jun 21, 2009 at 5:52 PM, <josef.pktd@gmail.com> wrote:
> On Sun, Jun 21, 2009 at 5:28 PM, <josef.pktd@gmail.com> wrote:
>> On Sun, Jun 21, 2009 at 4:13 PM, Celvin<read.beyond.data@gmx.net> wrote:
>>> josef.pktd@gmail.com @ Jun, 21, 2009 09:23PM
>>>
>>>> Did you try
>>>> UnivariateSpline.get_coeffs()
>>>
>>>> UnivariateSpline.get_knots()
>>>
>>> >From my experiments with interpolate splines, I would think this
>>>> provides what you want. But the documentation is still a bit sparse.
>>> Yes, I did try using UnivariateSpline. Apart from being "more off"
>>> than using splrep with a smoothing factor of 0, get_coeffs() also
>>> returns an 1-d array, with far too few coefficients.
>>>
>>> Assuming a signal with about 390 data points, I also expect about
>>> 390 coefficients (which is consistent with what I get using
>>> signal.cspline1d for example, but I need 4 arrays, not just one).
>>>
>>> I would expect coefficients to be a list of shape
>>>
>>> coeff = [ [a0, b0, c0, d0],
>>> [a1, b1, c1, d1],
>>> [a2, b1, c2, d2],
>>> ...
>>> ...
>>> [an, bn, cn, dn],
>>> ]
>>>
>>> but I am only able to obtain 1-d array of coefficients, no matter what
>>> function or module I use.
>>>
>>> http://www.physics.utah.edu/~detar/phys6720/handouts/cubic_spline/cubic_spline/node1.html
>>>
>>> ...is basically what I'm looking for. I used to do the matrix
>>> calculations for Si(x) myself using C++, but I'd prefer to somehow get the
>>> coefficients using numpy/scipy and not use an extension.
>>>
>>> Any further ideas?
>>
>> interpolate.spalde gives you the derivatives of as a 1 by 4 array
>> for each point
>> given it's a cubic function, it should be possible to recover the coefficients
>>
>> something like
>>
>> coeff = np.vstack(interpolate.spalde(x,tck))*[1,1,0.5,1/3.]
>
>
> I don't think that makes any sense. (sleep deprivation)
>
> Josef
>
>>
>> I have no idea yet, whether this actually does what you want.
>> If you find out how it works, it would be good to get an example since
>> figuring out the splines with the current documentation is pretty
>> difficult.
>
Here is a version that seems to work up to a precision of 1e-10. I
just wanted to see how it works, and didn't try for numerical
precision. I'm not sure it's useful, since the spline code calculates
the interpolation faster and with higher precision.
deriv2coeff(x2, tck):
uses the derivatives of the spline at an array of points to recover
the polynomial coefficients by direct calculation
Josef
"""
try_spline_coeff.py
"""
import numpy as np
from scipy import interpolate
def cubic(x,a):
return np.sum(a*np.c_[np.ones(x.shape), x, x**2, x**3],axis=-1)
def cubice(x,a):
return (a[:,0]+a[:,1]*x+a[:,2]*x**2+a[:,3]*x**3)
npoints = 51
x = np.linspace(0, 20, npoints)
y = np.sqrt(x) + np.sin(x) + 0.2* np.random.randn(npoints)
tck = interpolate.splrep(x, y)
x2 = np.linspace(0, 20, 200)
y2 = interpolate.splev(x2, tck)
def deriv2coeff(x2, tck):
'''get polynomial coefficients for spline at points x2'''
coeffr = np.vstack(interpolate.spalde(x2,tck))
coeffr[:,3] /= 6.
coeffr[:,2] = (coeffr[:,2] - 6*coeffr[:,3]*x2)/2.
coeffr[:,1] = coeffr[:,1] - 2*coeffr[:,2]*x2 - 3*coeffr[:,3]*x2**2
coeffr[:,0] = coeffr[:,0] - coeffr[:,1]*x2 - coeffr[:,2]*x2**2 -
coeffr[:,3]*x2**3
return coeffr
coeffr = deriv2coeff(x2, tck)
print cubic(x2[:10],coeffr[:10,:])
print cubice(x2[:10],coeffr[:10,:])
print y2[:10]
y2p = cubice(x2,coeffr)
print np.max(np.abs(y2-y2p))
y3s = interpolate.splev(x2+0.01, tck)
y3p = cubice(x2+0.01,coeffr)
print np.mean(np.abs(y3s-y3p)>1e-10)
y3s = interpolate.splev(x2-0.01, tck)
y3p = cubice(x2-0.01,coeffr)
print np.mean(np.abs(y3s-y3p)>1e-10)
x3 = tck[0][3:-3] # check at knots + 0.02
y3s = interpolate.splev(x3+0.02, tck)
coeffr0 = deriv2coeff(x3, tck)
y3p = cubice(x3+0.02,coeffr0)
print np.mean(np.abs(y3s-y3p)>1e-10)
More information about the SciPy-user
mailing list