如何显示UnivariateSpline生成的分段多项式显式表达式？

问题背景 #

我使用以下代码求解分段多项式函数：

import numpy as np
import matplotlib.pyplot as plt
import scipy
from scipy.interpolate import UnivariateSpline, splrep

x = np.array([0., 0.75, 1.8, 2.25, 3.75, 4.5, 6.45, 6.75, 
          7.5, 8.325, 10.875, 11.25, 12.525, 12.75, 15., 
          20.85, 21.])
y = np.array([2.83811035, 2.81541896, 3.14311655, 
          3.22373554, 3.43033456, 3.50433385, 
          3.66794514, 3.462296,   3.59480959, 
          3.56250726, 3.6209845,  3.63034523, 
          3.68238915, 3.69096892, 3.75560395, 
          3.83545191, 3.90419498])
          

k = 3  # 多项式次数
spl = UnivariateSpline(x, y, k=3, s=0.09)
xs = np.linspace(x.min(), x.max(), 100)
plt.plot(x, y, 'ro', ms=5)
plt.plot(xs, spl(xs), 'cyan', lw=5, alpha=0.3)

# 获取样条系数和节点
tck = (spl._data[8], spl._data[9], k)  # tck = (节点, 系数, 次数)
p = scipy.interpolate.PPoly.from_spline(tck)

# 绘制每段曲线并输出节点和系数
for idx, i in enumerate(range(k, len(spl.get_knots()) + k - 1)):
    xs = np.linspace(p.x[i], p.x[i + 1], 100)
    plt.plot(xs, np.polyval(p.c[:, i], xs - p.x[i]))
    print("knot ", p.x[i], " to ", p.x[i + 1])
    print("coeffs ", p.c[:, i], "\n")
    f0 = lambda x: p.c[0, i] * (x - p.x[i]) ** 3 + p.c[1, i] * (x - p.x[i]) ** 2 + p.c[2, i] * (x - p.x[i]) + p.c[3, i]
    f0 = lambda x: [p.c[:, i] * (x - p.x[i]) ** (3 - i) for i in range(k + 1)]
    print(f0)
plt.show()

遇到两个问题：

1. 执行print(f0)时，输出的是函数对象（如<function fit_spline1.<locals>.<lambda> at 0x0000028697B94F70>），希望显示具体的多项式表达式而非函数对象。
2. 参考一段代码后，运行得到16个基函数，但预期仅3个，不知道如何适配自己的场景：

from sympy import lambdify, bspline_basis_set
from sympy.abc import u

basis = bspline_basis_set(tck[2], tck[0],  u)
for i, b in enumerate(basis):
    print(f"Basis {i} :", b)

解决方案 #

1. 输出具体多项式表达式 #

直接定义lambda只能得到函数对象，要生成可读的多项式表达式，可以借助sympy库生成符号表达式，修改循环内的代码如下：

# 替换原循环中f0相关代码
import sympy as sp
x_sym = sp.symbols('x')

for idx, i in enumerate(range(k, len(p.x) - 1)):  # 修正循环范围，避免索引越界
    x0 = p.x[i]
    x1 = p.x[i+1]
    xs = np.linspace(x0, x1, 100)
    plt.plot(xs, np.polyval(p.c[:, i], xs - x0))
    
    # 输出区间和系数
    print(f"区间 [{x0:.3f}, {x1:.3f}]")
    print(f"系数（三次项→常数项）: {p.c[:, i]}\n")
    
    # 构建符号多项式表达式
    coeffs = p.c[:, i]
    poly_expr = coeffs[0]*(x_sym - x0)**3 + coeffs[1]*(x_sym - x0)**2 + coeffs[2]*(x_sym - x0) + coeffs[3]
    # 简化表达式
    simplified_expr = sp.simplify(poly_expr)
    print(f"多项式表达式：{simplified_expr}\n")

2. 理解B样条基函数数量问题 #

你得到16个基函数是正常的，B样条基函数数量由节点向量长度和样条次数决定，公式为：
基函数数量 = 节点向量长度 - 样条次数 - 1

你的样条是3次，tck[0]即节点向量的长度为20（包含两端重复的节点），因此20-3-1=16，和输出结果一致。

你需要的是每个分段区间的多项式表达式，而非B样条的基函数，因此无需使用bspline_basis_set，直接基于PPoly的结果生成即可，也就是上面第一部分的解决方案。

内容的提问来源于stack exchange，提问作者Natasha