1

假设我在 python 中有两个数组，我希望获得（并实际使用）这些点之间的三次样条插值。（即：我希望集成该功能）。我非常喜欢使用 numpy scipy 的方法。

我知道scipy.interpolate.interp1d。但是，这只允许我评估要点，例如非常简单的功能：

现在我可以做一些简单的事情： 

import numpy as np
import scipy.interpolate
import matplotlib.pyplot as plt

y = np.array([0,2,3,4,8,10,12,12,12,10,9,8,7,1,0,0,1,2])
x = np.array(range(len(y)))
xvals = np.linspace(0, len(y)-1, len(y)*100, endpoint = False)
func = scipy.interpolate.interp1d(x, y, kind = "cubic")
yvals = func(xvals)
plt.plot(xvals,yvals)
plt.plot(x,y, "o")

但是我希望进一步处理这个三次样条（即我需要得到积分）。对于手动操作我需要得到因素，所以：

a_i * x^3 + b_i * x^2 + c_i * x + d_i where i goes from 0 to n/3

（n = 元素的数量 - 这只是第 i 个三次的定义）

因此，我期望描述所有样条的元组（或二维数组）列表。- 或者一种获得第 i 个三次方的方法，并且真的，真的很想得到一个方便的“x-to-i”来找到我目前在哪个样条线。

（当然，后一个问题是在排序列表中简单地搜索大于引用的第一个值——如果需要，我可以很容易地手动完成）。

4

3 回答 3

1

只是一个可能会有所帮助的想法。 从文档中，您可以通过不同的方式获得三次样条插值，这可能对您有所帮助：

import numpy as np
import scipy.interpolate
import matplotlib.pyplot as plt

y = np.array([0,2,3,4,8,10,12,12,12,10,9,8,7,1,0,0,1,2])
x = np.array(range(len(y)))
xvals = np.linspace(0, len(y)-1, len(y)*100, endpoint = False)
func = scipy.interpolate.splrep(x, y, s=0)
yvals = scipy.interpolate.splev(xvals, func, der=0)

# display original vs cubic spline representation for security...
plt.figure()
plt.plot(x, y, 'x', xvals, yvals, x, y, 'b')
plt.legend(['Linear', 'Cubic Spline'])
plt.axis([-0.05, 20, -2, 20])
plt.title('Cubic-spline interpolation')
plt.show()

这使您可以通过以下方式访问系数

pp = scipy.interpolate.spltopp(func[0][1:-1],func[1],func[2])

#Print the coefficient arrays, one for cubed terms, one for squared etc
print(pp.coeffs)

它还在页面上提供了一个示例，说明如何使用此三次样条表示进行积分（我希望更改常数以适应您的情况 - 您的里程可能会有所不同）：

def integ(x, tck, constant=0):
    x = np.atleast_1d(x)
    out = np.zeros(x.shape, dtype=x.dtype)
    for n in xrange(len(out)):
        out[n] = scipy.interpolate.splint(0, x[n], tck)
    out += constant
    return out

const_of_integration = 0
yint = integ(xvals, func, const_of_integration)
于 2015-04-14T10:55:30.380 回答
1

对于插值，您可以使用scipy.interpolate.UnivariateSpline(..., s=0).

除其他外，它还具有集成方法。

编辑：s=0UnivariateSpline 构造函数的参数强制样条通过所有数据点。结果在 B 样条基础上，您可以使用get_coefs()和get_knots()方法获得节点和系数。格式与 netlib 上的 FITPACK、dierckx 中使用的格式相同。interp1d请注意， （依赖于splmakeATM）和UnivariateSpline（或者splrep/ ）内部使用的 tck 格式splev不一致。

EDIT2：您可以使用PPoly.from_spline获得样条曲线的分段多项式表示——但这与 iterp1d不一致。使用 splrep(..., s=0) 获得插值样条，然后转换结果。

于 2015-04-05T16:28:20.520 回答
0

该函数scipy.signal.cspline1d返回具有镜像对称边界条件的 b 样条项的系数。如果您需要在系数中执行额外的过滤，这很有用

于 2020-04-12T22:28:12.867 回答