使用 SymPy 通过给定点构造符号插值样条

Question

假设我从 R2 上定义的一些简单数据集开始：

DataPointsDomain = [0,1,2,3,4,5]
DataPointsRange =  [3,6,5,7,9,1]

使用 scipy 我可以使用以下方法制作惰性多项式样条：

ScipySplineObject = scipy.interpolate.InterpolatedUnivariateSpline( 
    DataPointsDomain, 
    DataPointsRange, 
    k = 1, )

sympy 中的等价对象是什么？？

SympySplineObject = ...???

（我想定义此对象并进行分析性的 sympy 操作，例如对 sympy 对象进行积分、导数等...）

Answer 1

在SymPy 1.1.1以上的版本中，包括current development version，有一个内置方法interpolating_spline，它有四个参数：样条阶数，变量，定义域值和范围值.

from sympy import *
DataPointsDomain = [0,1,2,3,4,5]
DataPointsRange =  [3,6,5,7,9,1]
x = symbols('x')
s = interpolating_spline(3, x, DataPointsDomain, DataPointsRange)

这个returns

Piecewise((23*x**3/15 - 33*x**2/5 + 121*x/15 + 3, (x >= 0) & (x <= 2)),
   (-2*x**3/3 + 33*x**2/5 - 55*x/3 + 103/5, (x >= 2) & (x <= 3)), 
   (-28*x**3/15 + 87*x**2/5 - 761*x/15 + 53, (x >= 3) & (x <= 5)))

这是通过给定点的 "not a knot" 三次样条。

旧答案

插值样条可以用 SymPy 构造，但这需要一些努力。 bspline_basis_set returns 方法是给定 x 值的 B 样条曲线的基础，但随后由您来找到它们的系数。

首先，我们需要节点列表，它与 x 值列表（下面的 xv）不完全相同。端点 xv[0] 和 xv[-1] 将出现 deg+1 次，其中 deg 是样条的度数，因为在端点处所有系数都会更改值（从某值变为零）。此外，一些接近它们的 x 值可能根本不会出现，因为那里的系数不会发生变化（"not a knot" 条件）。最后，对于偶数度样条（糟糕），内部结点位于数据点之间的中间。所以我们需要这个辅助函数：

from sympy import *
def knots(xv, deg):
    if deg % 2 == 1:
        j = (deg+1) // 2
        interior_knots = xv[j:-j]
    else:
        j = deg // 2
        interior_knots = [Rational(a+b, 2) for a, b in zip(xv[j:-j-1], xv[j+1:-j])]
    return [xv[0]] * (deg+1) + interior_knots + [xv[-1]] * (deg+1)

从 bspline_basis_set 方法获得 b 样条后，必须插入 x 值并形成一个线性系统，从中找到系数 coeff。最后，构造样条：

xv = [0, 1, 2, 3, 4, 5]
yv =  [3, 6, 5, 7, 9, 1]
deg = 3
x = Symbol("x")
basis = bspline_basis_set(deg, knots(xv, deg), x)
A = [[b.subs(x, v) for b in basis] for v in xv]
coeff = linsolve((Matrix(A), Matrix(yv)), symbols('c0:{}'.format(len(xv))))
spline = sum([c*b for c, b in zip(list(coeff)[0], basis)])
print(spline)

这个样条是一个 SymPy 对象。这是第 3 级：

3*Piecewise((-x**3/8 + 3*x**2/4 - 3*x/2 + 1, (x >= 0) & (x <= 2)), (0, True)) + Piecewise((x**3/8 - 9*x**2/8 + 27*x/8 - 27/8, (x >= 3) & (x <= 5)), (0, True)) + 377*Piecewise((19*x**3/72 - 5*x**2/4 + 3*x/2, (x >= 0) & (x <= 2)), (-x**3/9 + x**2 - 3*x + 3, (x >= 2) & (x <= 3)), (0, True))/45 + 547*Piecewise((x**3/9 - 2*x**2/3 + 4*x/3 - 8/9, (x >= 2) & (x <= 3)), (-19*x**3/72 + 65*x**2/24 - 211*x/24 + 665/72, (x >= 3) & (x <= 5)), (0, True))/45 + 346*Piecewise((x**3/30, (x >= 0) & (x <= 2)), (-11*x**3/45 + 5*x**2/3 - 10*x/3 + 20/9, (x >= 2) & (x <= 3)), (31*x**3/180 - 25*x**2/12 + 95*x/12 - 325/36, (x >= 3) & (x <= 5)), (0, True))/45 + 146*Piecewise((-31*x**3/180 + x**2/2, (x >= 0) & (x <= 2)), (11*x**3/45 - 2*x**2 + 5*x - 10/3, (x >= 2) & (x <= 3)), (-x**3/30 + x**2/2 - 5*x/2 + 25/6, (x >= 3) & (x <= 5)), (0, True))/45

你可以区分它，

spline.diff(x)

你可以集成它：

integrate(spline, (x, 0, 5))   #  197/3

你可以绘制它并看到它确实插入了给定的值：

plot(spline, (x, 0, 5))

我什至将它们一起绘制为 1、2、3 度：

免责声明：

上面给出的代码适用于 development version of SymPy，并且应该适用于 1.1.2+；以前版本的B样条方法有一个错误。
其中一些需要花费大量时间，因为 Piecewise 对象很慢。根据我的经验，基础建设需要最长的时间。

使用 SymPy 通过给定点构造符号插值样条

Construct a symbolic interpolating spline through given points using SymPy

python

sympy

curve-fitting

spline

旧答案