使用辛普森法计算 B 样条曲线弧长

1. 弧长的数学定义

设一条参数化的 B 样条曲线为 C(u)，其参数 u 的有效区间为 [a, b]。
该曲线的弧长 L 定义为：

L ＝ ∫ₐᵇ ‖C′(u)‖ du

其中：

• C′(u) 是曲线对参数 u 的一阶导数（即切向量），
• ‖C′(u)‖ 表示该切向量的欧几里得范数（即速度大小），
• a 和 b 是积分上下限（通常取 B 样条非零支撑区间的起止节点）。

由于该积分一般没有解析解，需采用数值积分方法求解。

2. 辛普森法基本公式（复合形式）

将区间 [a, b] 等分为 N 段（N 必须为偶数），步长为：

h ＝ (b − a) / N

记分点为：
u₀ ＝ a,
u₁ ＝ a ＋ h,
u₂ ＝ a ＋ 2h,
…,
uₙ ＝ b （其中 n ＝ N）

则弧长 L 可用复合辛普森公式近似为：

L ≈ (h / 3) ⋅ [ f(u₀) ＋ 4⋅f(u₁) ＋ 2⋅f(u₂) ＋ 4⋅f(u₃) ＋ 2⋅f(u₄) ＋ … ＋ 4⋅f(uₙ₋₁) ＋ f(uₙ) ]

其中函数 f(u) 定义为：

f(u) ＝ ‖C′(u)‖

更紧凑地写为：

L ≈ (h / 3) ⋅ [ f(u₀) ＋ f(uₙ) ＋ 4⋅Σₖ₌₁,₃,₅,…,ₙ₋₁ f(uₖ) ＋ 2⋅Σₖ₌₂,₄,₆,…,ₙ₋₂ f(uₖ) ]

3. 计算步骤

1. 确定参数区间 [a, b]：
   对于 p 次 B 样条，通常取第一个非零基函数起点到最后一个非零基函数终点，例如 [uₚ, uₘ₋ₚ]（m 为节点数）。

2. 选择偶数分割数 N（如 N ＝ 10, 20, 50, 100；越大精度越高）。

3. 计算步长：
   h ＝ (b − a) / N

4. 对每个分点 uᵢ ＝ a ＋ i⋅h（i ＝ 0, 1, …, N）：

   • 计算曲线在 uᵢ 处的一阶导数 C′(uᵢ)
   • 计算 f(uᵢ) ＝ ‖C′(uᵢ)‖

5. 按辛普森权重求和：

   • 首尾点（i ＝ 0 和 i ＝ N）权重为 1
   • 奇数索引点（i 为奇数）权重为 4
   • 偶数索引点（i 为偶数且 0 < i < N）权重为 2

6. 计算最终弧长：
   L ≈ (h / 3) × 总加权和

4. 示例（N ＝ 4）

若 N ＝ 4，则 h ＝ (b − a)/4，分点为 u₀, u₁, u₂, u₃, u₄。

弧长近似为：

L ≈ (h / 3) ⋅ [ f(u₀) ＋ 4f(u₁) ＋ 2f(u₂) ＋ 4f(u₃) ＋ f(u₄) ]

---

5. 注意事项

• N 必须为偶数，否则辛普森法不适用。
• 若曲线曲率变化剧烈（如高曲率拐点），建议局部加密分段或改用自适应辛普森法。
• 相比矩形法或梯形法，辛普森法具有 O(h⁴) 的收敛阶，精度更高。
• 对于光滑的 B 样条曲线（C² 连续），辛普森法通常只需 N ＝ 20～50 即可获得良好精度。

如需更高效率，也可结合自适应辛普森法（Adaptive Simpson’s Rule），在误差大的子区间自动递归细分。