Octave 中的三次样条实现

Question

我的大胆主张是，在 interp1(..., "spline") 中实现的三次样条的 Octave 实现不同于 "natural cubic spline" 算法，例如 Wolfram's Mathworld。我自己编写了后者的实现并将其与 interp1(..., "spline") 函数的输出进行比较，结果如下：

我发现当我尝试用 4 个点进行相同的比较时，解决方案也不同，而且，Octave 解决方案与将单个三次多项式拟合到所有四个点相同（实际上并没有产生分段样条对于三个间隔）。

我还尝试深入了解 Octave 的样条实现，发现它太过晦涩，无法在 5 分钟内阅读和理解。

我知道在实施三次样条时可以选择几个边界条件选项（"natural" 与 "clamped"）。我的实现使用 "natural" 边界条件（其中两个外部点的二阶导数设置为零）。

如果Octave的三次样条确实和标准的三次样条不同，那到底是什么？

编辑：

此处绘制了上面比较图中显示的两个解的二阶差分：

首先，在 Octave 的情况下似乎只有两个三次多项式：一个适用于前两个区间，一个适用于后两个区间。其次，他们显然没有使用 "natural" 样条曲线，因为极值处的二阶导数不趋于零。

此外，我认为我在中间（即第 3 个）点实现的二阶差分是零只是一个巧合，而不是算法所要求的.对一组不同的点重复此测试将 confirm/refute 这个。

Answer 1

不同的结束条件解释了您的实施与 Octave 的实施之间的差异。 Octave 使用 not-a-knot 条件（取决于输入）

见help spline

为了解释您的观察结果：由于非结条件，三阶导数在第 2 和第 (n-1) 次中断处是连续的，这就是为什么 Octave 的二阶导数看起来更少 'breaks'，因为它是前两段和后两段的连续直线。如果你看一下三阶导数，你可以更清楚地看到效果——三阶导数仅在第三次中断（中间）处不连续

x = 1:5;
y = rand(1,5);
xx = linspace(1,5);
pp = interp1(x, y, 'spline', 'pp');
yy = ppval(pp, xx);
dyy = ppval(ppder(pp, 3), xx);
plot(xx, yy, xx, dyy);

pp 数据结构也像这样

pp =

  scalar structure containing the fields:

    form = pp
    breaks =

       1   2   3   4   5

    coefs =

       0.427823  -1.767499   1.994444   0.240388
       0.427823  -0.484030  -0.257085   0.895156
      -0.442232   0.799439   0.058324   0.581864
      -0.442232  -0.527258   0.330506   0.997395

    pieces =  4
    order =  4
    dim =  1
    orient = first

Octave 中的三次样条实现

Cubic spline implementation in Octave

octave

cubic-spline