python - lagrange - coeficientes de interpolación spline en scipy

spline python

Question

Quiero calcular los coeficientes de una interpolación spline por scipy. En MATLAB:

x=[0:3];
y=[0,1,4,0];
spl=spline(x,y);
disp(spl.coefs);

y volverá:

ans =
-1.5000 5.5000 -3.0000 0
-1.5000 1.0000 3.5000 1.0000
-1.5000 -3.5000 1.0000 4.0000

Pero no puedo hacer eso por interpolate.splrep en scipy. ¿Puedes decirme cómo calcificarlo?

Answer 1

No estoy seguro de que haya alguna manera de obtener exactamente esos coeficientes de scipy. Lo que scipy.interpolate.splrep te da son los coeficientes para los nudos para una b-spline. Lo que le da spline de Matlab parece ser el coeficiente polinomial parcial que describe las ecuaciones cúbicas que conectan los puntos por los que pasa, lo que me lleva a pensar que la spline de Matlab es una spline basada en puntos de control como Hermite o Catmull-Rom en lugar de b-spline.

Sin embargo, scipy.interpolate.interpolate.spltopp proporciona una forma de obtener los coeficientes polinomiales parciales de una b-spline. Desafortunadamente, no parece funcionar muy bien.

>>> import scipy.interpolate
>>> x = [0, 1, 2, 3]
>>> y = [0, 1, 4, 0]
>>> tck = scipy.interpolate.splrep(x, y)
>>> tck
Out:
(array([ 0., 0., 0., 0., 3., 3., 3., 3.]),
array([ 3.19142761e-16, -3.00000000e+00, 1.05000000e+01,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00]),
3)
>>> pp = scipy.interpolate.interpolate.spltopp(tck[0][1:-1], tck[1], tck[2])
>>> pp.coeffs.T
Out:
array([[ -4.54540394e-322, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000],
[ -4.54540394e-322, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000],
[ -4.54540394e-322, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000],
[ 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000],
[ 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000]])

Tenga en cuenta que hay un conjunto de coeficientes por nudo, no uno para cada uno de los puntos originales pasados. Además, multiplicar los coeficientes por la matriz de base b-spline no parece ser muy útil.

>>> bsbm = array([[-1, 3, -3, 1], [ 3, -6, 3, 0], [-3, 0, 3, 0],
[ 1, 4, 1, 0]]) * 1.0/6
Out:
array([[-0.16666667, 0.5 , -0.5 , 0.16666667],
[ 0.5 , -1. , 0.5 , 0. ],
[-0.5 , 0. , 0.5 , 0. ],
[ 0.16666667, 0.66666667, 0.16666667, 0. ]])
>>> dot(pp.coeffs.T, bsbm)
Out:
array([[ 7.41098469e-323, -2.27270197e-322, 2.27270197e-322,
-7.41098469e-323],
[ 7.41098469e-323, -2.27270197e-322, 2.27270197e-322,
-7.41098469e-323],
[ 7.41098469e-323, -2.27270197e-322, 2.27270197e-322,
-7.41098469e-323],
[ 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000],
[ 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000]])

El Paquete de Polinomios Piecewise de FORTRAN, PPPack , tiene un comando bsplpp que se convierte de B-spline a forma de polinomio por partes, que puede servir a sus necesidades. Desafortunadamente, no hay un contenedor de Python para PPPack en este momento.

Answer 2

Los documentos en scipy.interpolate.splrep dicen que puedes obtener los coeficientes:

Returns:
tck : tuple
(t,c,k) a tuple containing the vector of knots, the B-spline coefficients, and the degree of the spline.

Answer 3

Aquí es cómo podría obtener resultados similares a MATLAB:

>>> from scipy.interpolate import PPoly, splrep
>>> x = [0, 1, 2, 3]
>>> y = [0, 1, 4, 0]
>>> tck = splrep(x, y)
>>> tck
Out: (array([ 0., 0., 0., 0., 3., 3., 3., 3.]),
array([ 3.19142761e-16, -3.00000000e+00, 1.05000000e+01,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00]),
3)
>>> pp = PPoly.from_spline(tck)
>>> pp.c.T
Out: array([[ -1.50000000e+00, 5.50000000e+00, -3.00000000e+00,
3.19142761e-16],
[ -1.50000000e+00, 5.50000000e+00, -3.00000000e+00,
3.19142761e-16],
[ -1.50000000e+00, 5.50000000e+00, -3.00000000e+00,
3.19142761e-16],
[ -1.50000000e+00, 5.50000000e+00, -3.00000000e+00,
3.19142761e-16],
[ -1.50000000e+00, -8.00000000e+00, -1.05000000e+01,
0.00000000e+00],
[ -1.50000000e+00, -8.00000000e+00, -1.05000000e+01,
0.00000000e+00],
[ -1.50000000e+00, -8.00000000e+00, -1.05000000e+01,
0.00000000e+00]])