logo标签列表

Scipy中样条插值的系数

pythonmatlabscipy
12

我想使用scipy计算样条插值的系数。 在MATLAB中：

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

它会返回：

ans =

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

但是我无法通过Scipy中的interpolate.splrep进行计算。你能告诉我如何计算吗？

1
你可以轻松地获得节点系数。这样可以吗？ - brentlance
它不像MATLAB的结果那样返回系数。 - Gochit
4个回答
8
我不确定是否有办法从scipy获得确切的系数。scipy.interpolate.splrep给出的是B样条曲线节点的系数。而Matlab的spline似乎给出了描述连接传入点的三次方程的偏导多项式系数。这让我相信Matlab的spline是一种基于控制点的样条，如Hermite或Catmull-Rom，而非B样条。
然而，scipy.interpolate.interpolate.spltopp提供了一种获取B样条的偏导多项式系数的方法。不幸的是，它似乎并不起作用。
>>> 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]])

请注意，每个节点只有一个系数集，而不是针对传入的每个原始点都有一个系数。此外，将系数乘以B样条基矩阵似乎并不是很有帮助。
>>> 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]])

FORTRAN分段多项式包PPPack提供了一个命令bsplpp，可以将B样条转换为分段多项式形式，也许可以满足您的需求。不幸的是，目前还没有Python的PPPack包装器。

谢谢，你的回答很有帮助，增加了我的知识。但是我通过另一个样条插值的Python脚本解决了这个问题。 - Gochit
没问题。很高兴能帮到你。你用了哪个脚本来解决这个问题？ - brentlance
1
我在《Python工程数值方法》中使用了样条插值脚本。 - Gochit
我不知道这本书的存在。谢谢！ - brentlance
1
还有from_spline，pp.from_spline(tck).c 似乎具有更合理的系数。 - dashesy
3

如果您安装了Scipy版本>= 0.18.0，您可以使用来自Scipy插值模块的CubicSpline函数进行三次样条插值。

您可以通过在Python中运行以下命令来检查Scipy版本:

#!/usr/bin/env python3
import scipy
scipy.version.version

如果您的scipy版本为 0.18.0 或更高，则可以运行以下代码示例进行三次样条插值：
#!/usr/bin/env python3

import numpy as np
from scipy.interpolate import CubicSpline

# calculate 5 natural cubic spline polynomials for 6 points
# (x,y) = (0,12) (1,14) (2,22) (3,39) (4,58) (5,77)
x = np.array([0, 1, 2, 3, 4, 5])
y = np.array([12,14,22,39,58,77])

# calculate natural cubic spline polynomials
cs = CubicSpline(x,y,bc_type='natural')

# show values of interpolation function at x=1.25
print('S(1.25) = ', cs(1.25))

## Aditional - find polynomial coefficients for different x regions

# if you want to print polynomial coefficients in form
# S0(0<=x<=1) = a0 + b0(x-x0) + c0(x-x0)^2 + d0(x-x0)^3
# S1(1< x<=2) = a1 + b1(x-x1) + c1(x-x1)^2 + d1(x-x1)^3
# ...
# S4(4< x<=5) = a4 + b4(x-x4) + c5(x-x4)^2 + d5(x-x4)^3
# x0 = 0; x1 = 1; x4 = 4; (start of x region interval)

# show values of a0, b0, c0, d0, a1, b1, c1, d1 ...
cs.c

# Polynomial coefficients for 0 <= x <= 1
a0 = cs.c.item(3,0)
b0 = cs.c.item(2,0)
c0 = cs.c.item(1,0)
d0 = cs.c.item(0,0)

# Polynomial coefficients for 1 < x <= 2
a1 = cs.c.item(3,1)
b1 = cs.c.item(2,1)
c1 = cs.c.item(1,1)
d1 = cs.c.item(0,1)

# ...

# Polynomial coefficients for 4 < x <= 5
a4 = cs.c.item(3,4)
b4 = cs.c.item(2,4)
c4 = cs.c.item(1,4)
d4 = cs.c.item(0,4)

# Print polynomial equations for different x regions
print('S0(0<=x<=1) = ', a0, ' + ', b0, '(x-0) + ', c0, '(x-0)^2  + ', d0, '(x-0)^3')
print('S1(1< x<=2) = ', a1, ' + ', b1, '(x-1) + ', c1, '(x-1)^2  + ', d1, '(x-1)^3')
print('...')
print('S5(4< x<=5) = ', a4, ' + ', b4, '(x-4) + ', c4, '(x-4)^2  + ', d4, '(x-4)^3')

# So we can calculate S(1.25) by using equation S1(1< x<=2)
print('S(1.25) = ', a1 + b1*0.25 + c1*(0.25**2) + d1*(0.25**3))

# Cubic spline interpolation calculus example
    #  https://www.youtube.com/watch?v=gT7F3TWihvk
2

以下是我如何获得类似于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]])
你如何使用coeff.形成方程？你能帮忙解释为什么数组中的coeff.是重复的吗？ - Legolas
@Legolas，你可以使用sympy，例如from sympy import poly。你应该创建一个新的问题，并更清楚地表达你想要什么。 - dashesy
1

在scipy.interpolate.splrep的文档中提到，您可以获取系数：

Returns:

  tck : tuple

  (t,c,k) a tuple containing the vector of knots, the B-spline coefficients, and the degree of the spline.
我尝试了：tck=interpolate.splrep(x,y,k=3,s=0) 它返回系数：array([ -4.09317213e-17, -3.00000000e+00,1.05000000e+01,0.00000000e+00,0.00000000e+00, 0.00000000e+00,0.00000000e+00,0.00000000e+00])。 这与我的MATLAB结果不同。 - Gochit
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接