de Boor 样条求值 | Codeboy's Blog

形式陈述

输入为次数 $p$ 、夹持非降节点向量 $t_0,\ldots,t_{N+p}$ 、系数 $c_0,\ldots,c_{N-1}$ ，输出

s(x)=\sum_{i=0}^{N-1}c_iB_{i,p}(x).

先检查两端各重 $p+1$ 次、内部重数不超过 $p+1$ 、 $a<b$ 、系数数量正确以及数值有限。本文只接受 $x\in[a,b]$ ；区间外不自动外推。

若 $x<b$ ，找唯一非零跨度 $t_k\le x<t_{k+1}$ ，即相等节点块之后的那一段；若 $x=b$ ，令 $k=N-1$ ，按左极限处理。复制 $d_i^{(0)}=c_i$ ， $i=k-p,\ldots,k$ 。对 $r=1,\ldots,p$ ，计算

\begin{aligned} \alpha_{i,r}&=\frac{x-t_i}{t_{i+p-r+1}-t_i},\\ d_i^{(r)}&=(1-\alpha_{i,r})d_{i-1}^{(r-1)} +\alpha_{i,r}d_i^{(r-1)},\qquad i=k-p+r,\ldots,k. \end{aligned}

返回 $d_k^{(p)}$ 。若在同一数组覆盖，每轮必须按 $i=k,k-1,\ldots,k-p+r$ 从右向左更新，避免读到本轮的新值。

有效跨度保证每个实际使用的分母为正： $i\le k$ 而 $i+p-r+1\ge k+1$ ，因此它至少跨越 $t_k<t_{k+1}$ 。这里不需要把有效步骤的 $0/0$ 人为改成零；若出现零分母，应检查跨度或下标是否有误。基函数递推里省略零支撑项的约定，不能不加区分地套到系数三角表。

直觉

Cox–de Boor 递推从低阶基构造高阶基。将 $s$ 展开后反向收集同一低阶基的系数，就得到上面的 $d_i^{(1)}$ ；再收集一次得 $d_i^{(2)}$ 。最后只剩跨度 $[t_k,t_{k+1})$ 的零次基，其值为一，旁边的系数就是所求值。这给出逐层不变式

s(x)=\sum_{i=k-p+r}^{k}d_i^{(r)}B_{i,p-r}(x),

其中每层只在当前查询点使用该等式。第一层到下一层的系数合并恰是递推式，末层为单项。

在合法跨度上 $0\le\alpha_{i,r}\le1$ 。每一步只混合两个邻系数，解释了凸包性质与局部计算。一次完整合法性预检需 $O(N+p)$ 时间，可在固定节点和系数后缓存。预检通过的单次查询定位需 $O(\log N)$ ，之后 $p(p+1)/2$ 次混合，合计 $O(\log N+p^2)$ 时间、 $O(p)$ 额外空间；它不必先生成全部 $N$ 个基值。

例子与边界

取 $p=2$ 、 $t=(0,0,0,1,1,2,3,3,3)$ 、 $c=(0,1,2,0,1,3)$ 。在 $x=3/2$ ， $k=4$ ，初行是

(d_2^{(0)},d_3^{(0)},d_4^{(0)})=(2,0,1).

第一轮 $\alpha_{3,1}=1/2$ 、 $\alpha_{4,1}=1/4$ ，得到 $(d_3^{(1)},d_4^{(1)})=(1,1/4)$ ；第二轮 $\alpha_{4,2}=1/2$ ，得到 $d_4^{(2)}=5/8$ 。这与逐基求值的 $2(1/4)+0(5/8)+1(1/8)$ 相同。

在二重节点 $x=1$ ，仍取 $k=4$ 。第一轮两个权重都为零，得到 $(2,0)$ ；第二轮权重为零，返回二。这是合法的零权重，不是零分母。若误把 $[t_3,t_4]=[1,1]$ 当作有效跨度，才会制造无意义的除零。

在 $x=3$ ，选末跨度 $k=5$ ，相关权重都为一，返回 $c_5=3$ 。次数零时没有混合步骤，直接返回当前跨度的系数。内部满重数允许左右值不同：例如 $p=1,t=(0,0,1,1,2,2),c=(0,1,4,5)$ ，左极限为一，而算法在 $x=1$ 返回右侧值四。

推论与应用

向量系数可逐分量执行同一表，得到参数曲线上的点。只有一个跨度时，这个三角表就是 Bernstein 表示下的 de Casteljau 求值；局部线性混合是两者共有的计算结构。

浮点凸组合避免展开大幂系数，但结果接近零时，相对误差仍可能很大；节点差极小时，减法与除法也需要尺度检查。验收可比较同一输入的基递推、分段多项式和三角求值，并检查端点、节点左右两侧以及非法输入。通过这些样例说明实现符合所定接口，不等于已经证明任意浮点输入都没有误差。

参考资料

Carl de Boor，On calculating with B-splines，1972，§1 式 (16)–(18)、§2：局部系数三角求值。
Carl de Boor，B(asic)-Spline Basics，§9，Algorithm 9：从降低阶数导出求值算法。

形式陈述 ​

反向链接

形式陈述