C#版Cubic Spline Interpolation代码公式原理相关疑问

1 三对角矩阵构造逻辑解释 #

你看到的系数和右端项公式都来自三次样条的三弯矩方程推导，这是自然三次样条最常用的求解形式：

• 我们将每个插值节点i处的样条二阶导数记为M_i，作为待求解的变量
• 相邻节点间距记为h[i] = x[i] - x[i-1]，节点处的原始值记为_values[i]
• 利用相邻两段样条在节点i处的一阶导数连续条件，可以推导出标准三弯矩方程：
  h[i] * M[i-1] + 2*(h[i] + h[i+1]) * M[i] + h[i+1] * M[i+1] = 6 * [ (_values[i+1] - _values[i])/h[i+1] - (_values[i] - _values[i-1])/h[i] ]
• 代码实现时为了简化计算，直接把方程两边同时除以6，就得到了你看到的系数：
  • 主对角线项diag[i] = (h[i] + h[i+1])/3：对应方程除以6后M[i]的系数2*(h[i]+h[i+1])/6 = (h[i]+h[i+1])/3
  • 上对角线项sup[i] = h[i+1]/6：对应方程除以6后M[i+1]的系数
  • 次对角线项sub[i] = h[i]/6：对应方程除以6后M[i-1]的系数
  • 右端项_a[i]就是方程右边的原始项，物理意义是相邻两段样条的平均斜率差值，代码里先把右端项存在_a数组里，后续求解三对角矩阵后，_a数组会被覆盖为每个节点的二阶导数M[i]的解。

2 GetValue方法返回值逻辑解释 #

这部分是基于M参数的三次样条插值公式的展开形式：
首先定位到插值点key所在的区间是[x[gap-1], x[gap]]，区间长度h = _h[gap]，代码里的x1 = key - x[gap-1]，x2 = x[gap] - key，显然满足x1 + x2 = h。
标准的基于二阶导数M的样条插值公式为：

S(x) = M[gap-1] * x2^3/(6h) + M[gap] * x1^3/(6h) + (_values[gap-1] - M[gap-1] * h²/6)*x2/h + (_values[gap] - M[gap] * h²/6)*x1/h

把上述公式合并同类项、因式分解整理后，就会得到代码里的形式：

[ (-M[gap-1]/6 * (x2 + h) * x1 + _values[gap-1]) * x2 + (-M[gap]/6 * (x1 + h) * x2 + _values[gap]) * x1 ] / h

和你贴的代码完全一致：代码里的_a[gap-1]、_a[gap]就是之前求解得到的节点二阶导数M[gap-1]、M[gap]。 #

内容的提问来源于stack exchange，提问作者Siong Jie Ting