如何手动评估由scipy生成的三次B样条（基于系数与节点）

make_interp_spline默认生成三次B样条（阶数k=3），其系数是B样条基函数的权重，而非直接的多项式系数，因此不能直接用节点间相对距离代入多项式计算。以下是手动评估的核心逻辑和伪代码：

核心前提 #

• t：完整节点数组（包含开头/结尾的重复节点，比如你提供的例子中，原节点数组可能为[0,0,0,0,2,3,4,5,5,5,5]，三次样条要求首尾各有k+1个重复节点）
• c：系数数组（你提供的[0., 34.66324, 86.65809, 155.98457, 207.97942]）
• k：样条阶数，默认k=3

手动评估步骤（伪代码） #

1. 定位待评估点的区间 #

找到最大的索引i，满足t[i] ≤ x < t[i+1]。边界点x = t[-1]直接取最后一个区间（i = len(t)-2）。用二分查找效率更高，适合Numba编译：

def find_interval(x, t):
    if x >= t[-1]:
        return len(t) - 2
    low, high = 0, len(t) - 1
    while low < high:
        mid = (low + high) // 2
        if x < t[mid]:
            high = mid
        else:
            low = mid + 1
    return low - 1

2. 用De Boor算法计算B样条基函数值 #

三次样条需要计算当前区间对应的4个基函数值，通过递推公式实现（避免递归，适配Numba）：

def compute_cubic_basis(x, t, i):
    # 取当前区间对应的4个基函数的节点范围
    j_start = i - 3
    b = [0.0]*4

    # 零次基函数（k=0）
    for idx in range(4):
        tj = t[j_start + idx]
        tj1 = t[j_start + idx + 1]
        if tj <= x < tj1:
            b[idx] = 1.0
        elif x == t[-1] and (j_start + idx + 1) == len(t):
            b[idx] = 1.0

    # 递推一次基函数（k=1）
    new_b = [0.0]*3
    for idx in range(3):
        denom1 = t[j_start + idx + 1] - t[j_start + idx]
        term1 = ((x - t[j_start + idx])/denom1) * b[idx] if denom1 != 0 else 0.0
        denom2 = t[j_start + idx + 2] - t[j_start + idx + 1]
        term2 = ((t[j_start + idx + 2] - x)/denom2) * b[idx+1] if denom2 != 0 else 0.0
        new_b[idx] = term1 + term2
    b = new_b

    # 递推二次基函数（k=2）
    new_b = [0.0]*2
    for idx in range(2):
        denom1 = t[j_start + idx + 2] - t[j_start + idx]
        term1 = ((x - t[j_start + idx])/denom1) * b[idx] if denom1 != 0 else 0.0
        denom2 = t[j_start + idx + 3] - t[j_start + idx + 1]
        term2 = ((t[j_start + idx + 3] - x)/denom2) * b[idx+1] if denom2 != 0 else 0.0
        new_b[idx] = term1 + term2
    b = new_b

    # 递推三次基函数（k=3）
    new_b = [0.0]*1
    denom1 = t[j_start + 3] - t[j_start]
    term1 = ((x - t[j_start])/denom1) * b[0] if denom1 != 0 else 0.0
    denom2 = t[j_start + 4] - t[j_start + 1]
    term2 = ((t[j_start + 4] - x)/denom2) * b[1] if denom2 != 0 else 0.0
    new_b[0] = term1 + term2

    # 补全另外三个三次基函数值（重复递推逻辑，此处简化为完整计算）
    # 实际可合并为通用递推函数，此处为清晰拆分
    return [b0, b1, b2, b3]  # 对应B_{i-3,3}, B_{i-2,3}, B_{i-1,3}, B_{i,3}

3. 计算样条值 #

取对应区间的4个系数，与基函数值加权求和：

def evaluate_spline(x, t, c):
    i = find_interval(x, t)
    basis = compute_cubic_basis(x, t, i)
    # 取i-3到i的系数
    coeffs = c[i-3:i+1]
    return sum(coeff * b for coeff, b in zip(coeffs, basis))

Numba适配优化版 #

将上述逻辑合并为Numba可编译的函数，消除Python对象开销：

import numba as nb

@nb.njit
def numba_evaluate_spline(x, t, c):
    k = 3
    m = len(t)
    # 定位区间
    if x >= t[-1]:
        i = m - 2
    else:
        low, high = 0, m - 1
        while low < high:
            mid = (low + high) // 2
            if x < t[mid]:
                high = mid
            else:
                low = mid + 1
        i = low - 1
    
    idx_start = i - k
    # 初始化基函数
    b = nb.typed.List([0.0]*(k+1))
    for j in range(k+1):
        tj = t[idx_start + j]
        tj1 = t[idx_start + j + 1]
        if tj <= x < tj1:
            b[j] = 1.0
        elif x == t[-1] and (idx_start + j + 1) == m:
            b[j] = 1.0
    
    # 递推计算高次基函数
    for d in range(1, k+1):
        new_b = nb.typed.List([0.0]*(k+1 - d))
        for j in range(len(new_b)):
            tj = t[idx_start + j]
            tj_d = t[idx_start + j + d]
            denom1 = tj_d - tj
            term1 = ((x - tj)/denom1) * b[j] if denom1 != 0 else 0.0
            
            tj1 = t[idx_start + j + 1]
            tj1_d = t[idx_start + j + 1 + d]
            denom2 = tj1_d - tj1
            term2 = ((tj1_d - x)/denom2) * b[j+1] if denom2 != 0 else 0.0
            
            new_b[j] = term1 + term2
        b = new_b
    
    # 加权求和
    value = 0.0
    for j in range(len(b)):
        value += c[idx_start + j] * b[j]
    return value

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