65.9K
CodeProject 正在变化。 阅读更多。
Home

从多项式到自然样条的插值

Mosi_62
starIconstarIconstarIconstarIcon
emptyStarIcon
starIcon

4.96/5 (48投票s)

2015年3月5日

CPOL

12分钟阅读

viewsIcon

58136

downloadIcon

5872

从多项式到自然样条的插值

引言

插值的问题最近来到我这里。所以我开始研究这个课题。

我实现了多项式、拉格朗日、牛顿和自然样条算法，并开始将它们进行比较。令我非常惊讶的是，我发现花哨的样条算法在任何情况下都不是最佳解决方案。对于某些波形，即使是多项式插值也比样条插值效果更好。

但首先，我将对这些算法进行一些解释...

多项式插值

多项式插值算法为 n 个支撑点构建一个 n 次多项式，如下所示：

其中 x 和 y 是一个支撑点的坐标。对于 n 个支撑点，我们得到 n 个这样的方程，从 x1, y1 到 xn, yn。因此，算法基本上需要建立一个 n*n 的方程矩阵，并通过高斯算法求解。这将得到参数 a0 到 an。有了这些参数，就可以计算出任意 xp 对应的 yp 值。

多项式插值是这 4 种算法中最容易实现的。但它在精度方面很快就会达到极限。如果支撑点之间的 deltaX 太小或太大，即使有 10 个支撑点，高斯算法在处理方程矩阵的配置时也会遇到问题。通过缩放 deltaX 的值（增大或减小），使其接近 1，可以有助于求解矩阵方程。我在进行大部分插值时使用了 deltaX 从 0.1 到 0.4。这效果相当不错。但问题依然存在。当支撑点超过 10 个时，多项式插值就不再可靠了。如下图所示：

使用 f=1/x^2 和 20 个支撑点进行多项式插值尝试。该现象早在 16 个支撑点时就已经开始出现。

但这还不是全部。

在方程矩阵中，我们从以下方程开始：

继续到 y2、y3 等方程。这意味着在第一行，由于 x1，最上面的值最小，而下面的值越来越大，高斯算法不喜欢这样。将方程矩阵颠倒填充，例如...

...

...可以显著提高多项式算法的精度。这样，对于 deltaX 为 0.5 的情况，20 个支撑点就可以被正确地插值。

    // fill matrix upside down
    for (i = 0; i < order; i++)
    {
        m.y[i] = y_k[order - 1 - i];
        for (k = 0; k < order; k++)
        {
            m.a[i, k] = Math.Pow(x_k[order - 1 - i], k);
        }
    }
...

仅仅这样填充矩阵，而不是

    // fill matrix reular
    for (i = 0; i < order; i++)
    {
        m.y[i] = y_k[i];
        for (k = 0; k < order; k++)
        {
            m.a[i, k] = Math.Pow(x_k[i], k);
        }
    }
...

就可以显著提高精度。

拉格朗日插值

拉格朗日为每个支撑值（包括插值点 xp）构建一个多项式。对于 n 个支撑值，这就产生了 n 个多项式 L1..Ln。

xp 是插值变量，需要计算其对应的 yp 值。

yp 现在是所有 Lk* yk 项的总和。

在 Lk 的多项式中可以看出，如果 xp = xk，则 Lk 的多项式变为 1，yp 变为 yk。xp 越接近其中一个支撑点，该点获得的权重就越大。

对于一个插值点，插值过程在一个循环中完成，包括计算所有 Lk 元素并将 Lk * yk 相加。

    for (i = 0; i < order; i++)
    {
        nominator = 1.0;
        denominator = 1.0;
        /* calculate the Lk elements */
        for (k = 0; k < order; k++)
        {
            if (k != i)
            {
                /* nominator and denominator for one L element */
                nominator = nominator * (xp - x_k[k]);
                denominator = denominator * (x_k[i] - x_k[k]);
            }
        }
        /* put every thing thogether to yp*/
        if (denominator != 0.0)
            yp = yp + y_k[i] * nominator / denominator;
        else
            yp = 0.0;
    }var i = 0;
...

此循环计算一个 yp 插值值。

当支撑点数量超过 10 个时，拉格朗日算法比多项式算法更可靠。它在 deltaX 很小时也会遇到问题。但是，插值 60 个支撑点，deltaX = 0.1，都没有问题。

牛顿插值

牛顿使用另一个多项式来插值 yp。

他计算系数 b0 到 bn 的方法如下：

在这种形式下，计算过程并不太直观。但如果我们把 |XkXk-1|… 项放入一个表格（这里 n=4），情况就会变得更清楚。

 

我们首先需要计算 |XkXk-1| 项。然后计算 |XkXk-1Xk-2|… 项，以此类推。尽管我们只需要每一列最顶端的值来计算 b0 到 b4，但我们必须计算所有元素。这看起来相当复杂。但通过递归函数，它会变得相当简单。

    void CalcElements(double[] x, int order, int step)
    {
        int i;
        double[] xx;
        if (order >= 1)
        {
            xx = new double[order];
            for (i = 0; i < order-1; i++)
            {
                xx[i] = (x[i + 1] - x[i]) / (x_k[i + step] - x_k[i]);
            }
            b[step - 1] = x[0];
            CalcElements(xx, order - 1, step + 1);
        }
    }

这个函数通过包含支撑点的数组和作为阶数的支撑点数量来调用，并从第一步开始。

现在可以像这样计算一个 yp 插值值：

    for (i = 1; i < order; i++)
    {
        tempYp = b[i];
        for (k = 0; k < i; k++)
        {
            tempYp = tempYp * (xp - x_k[k]);
        }
        yp = yp + tempYp;
    }

牛顿算法对小的 deltaX 非常不敏感，并且在 60 个支撑点、deltaX=0.001 的情况下都能很好地工作。我认为它是这 3 种算法中最好的。

对于更复杂的问题，则需要使用自然样条插值。

自然样条

对于样条插值，为每个支撑点之间的区间计算一个插值函数。

对于这些区间中的每一个，都计算一个三次多项式，如下所示：

对于一个从 xi 开始到 xi+1 结束的区间，以 x 作为插值变量。

为了确保所有这些区间的波形能够平滑连接，需要满足一些条件：

对于 i = 1 到 n-1

对于 i = 2 到 n-1

对于 i = 2 到 n-1

对于 i = 2 到 n-1

从这些条件和替换 hi = xi+1 - xi，可以推导出参数 a, b, c 和 d。

从条件 a，我们得到 ai = yi。

从 b，我们得到：

c 得到：

这两个方程组合在一起，得到一个关于 ci 的 n-2 阶矩阵方程，其矩阵为：

 

和解向量：

通过高斯算法求解这个矩阵方程，得到 c2 到 cn-1。第一个和最后一个 c 都设置为 0。

现在只剩下 di，我们从条件 d 得到：

有了这些参数，就可以对支撑点之间的每个区间进行插值，最后将这些片段组合成一个图形。

自然样条算法在 100 个支撑点까지 都能正常工作。然后它在精度上会遇到问题。但它无疑是最复杂的算法。

样条速度改进

为了计算 h 参数，我们需要求解一个具有以下结构的矩阵方程：

 

x1

x2

x3

x4

 

a1

b1

 

 

= d1

c2

a2

b2

 

= d2

 

c3

a3

b3

= d3

 

 

c4

a4

= d4

n=4 时矩阵方程

这是一个三对角矩阵（至少在我们不希望插值周期函数时是如此）。这种特殊矩阵方程可以通过将方程 Ax = D 的 A 矩阵分解为左三角矩阵和右三角矩阵（A = LR）来比高斯消元法更快地求解。

a1

b1

 

 

 

1

 

 

 

 

m1

r1

 

 

c2

a2

b2

 

 

l1

1

 

 

 

 

m2

r2

 

 

c3

a3

b3

=

 

l2

1

 

*

 

 

m3

r3

 

 

c4

a4

 

 

 

l3

1

 

 

 

 

m4

 

比较参数可得：

 

a1 = m1

b1 = r1

c2 = l1 * m1

a2 = l1 * r1 + m1

b2 = r2

c3 = l2 * m2

a3 = l2 * r2 + m2

b3 = r3

c4 = l3 * m3

a4 = l3 * r3 + m3

b4 = r4

 

由此可得：

            li = ci+1 / mi

            mi+1 = ai+1 – li * bi

            对于 i = 1 到 n-1

 

为了求解这个方程，我们首先将方程 LRx = D 扩展为 y，并说 LRxy = Dy（y 不是最上面方程的解向量）。这可以分离为 Ly = D 和 Rx = y。

Ly = D 看起来像：

1

 

 

 

 

y1

 

d1

l1

1

 

 

 

y2

 

d2

 

l2

1

 

*

y3

=

d3

 

 

l3

1

 

y4

 

d4

 

求解 Ly = D 得到：

 

y1 = d1
  • y1 = d1

y1 * l1 + y2 = d2
  • y2 = d2 - y1 * l1

y2 * l2 + y3 = d3
  • y3 = d3 – y2 * l2

y3 * l3 + y4 = d4
  • y4 = d4 – y3 * l3

 

Rx = y 看起来像：

m1

r1

 

 

 

x1

 

y1

 

m2

r2

 

 

x2

 

y2

 

 

m3

r3

*

x3

=

y3

 

 

 

m4

 

x4

 

y4

 

求解 Rx = y 得到：

x1 * m1 + x2 * r1 = y1
  • x1 = (y1 – x2 * r1) / m1

x2 * m2 + x3* r2 = y2
  • x2 = (y2 – x3 * r2) / m2

x3 * m3 + x4 * r3 = y3
  • x3 = (y3 – x4 * r3) / m3

x4 * m4 = y4
  • x4 = y4 / m4

 

有了这些，我们基本上需要实现 3 个循环来求解三对角矩阵方程。

 

第一个循环

m1 = a1
对于 i = 1 到 n-1
            li = ci+1 / mi
            mi+1 = ai+1 – li * bi

 

第二个循环

y1 = d1
对于 i = 2 到 n
            yi = di – li-1 * yi-1

 

第三个循环

xn = yn / mn
对于 I = n-1 到 1
            xi = (yi - bi * xi+1)/mi

 

现在，我们可以将第一个和第二个循环合并为一个，如果我们分离第一个循环的第一项和第二个循环的最后一项。

 

m1 = a1
l1 = c2 / m1
m2 = a2 – l1 * b1

对于 i = 2 到 n-1
            li = ci+1 / mi
            mi+1 = ai+1 – li * bi
            yi = di – li-1 * yi-1

yn = dn – ln-1 * yn-1

 

因此，只剩下 2 个循环就可以求解整个矩阵方程，这比高斯算法所需的计算时间少一半多。

我在示例项目 Splien_LR 的 SolveLR() 函数中实现了这个算法 :-)

 

比较

 

这 4 种算法的比较有点令人惊讶。多项式、拉格朗日和牛顿算法在某些条件下容易出现过冲。这一点到处都有说明，而且确实如此。只要支撑点的值上下波动，插值波形就会波动得更厉害。如 6 个支撑点的示例形状所示：

是

0.5

11.7

14.8

4.0

2.2

0.2

x

0.1

0.2

0.3

0.4

0.5

0.6

所有 3 种算法都显示出一些过冲。

样条插值产生了更好的结果。

我在一本书中找到的形状是 4*exp-(x-0.4)^2。

是

0.0005

0.073

1.472

4.0

1.472

0.073

0.0005

x

0.1

0.2

0.3

0.4

0.5

0.6

0.7

这似乎对多项式、牛顿和拉格朗日算法来说实在太难了。

样条算法处理得更好。

样条曲线倾向于更直接地接近支撑点，并以更小的半径弯曲。这看起来相当不错。但它并非在所有情况下都是最佳解决方案。

像 1/x 这样的波形，带有以下点：

是

1.0

0.5

0.3333

0.25

0.2

0.1666

0.143

0.125

x

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

变化相当大。

在这里，多项式、牛顿和拉格朗日算法会产生一个平滑、圆润的形状。那么样条呢？

令人惊喜的是：这更像一个香蕉形状。

所以，让我们更进一步，看看 1/x^2。

是

1.0

0.25

0.111

0.0625

0.04

0.0277

0.0204

0.0156

x

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

多项式、牛顿和拉格朗日算法再次产生了平滑的形状。

以及样条曲线：

好的，它不再像香蕉了。:)

Using the Code

我尽可能地实现了所有 4 种算法。有一个全局常量“order”定义了支撑点的数量。这个数量可以设置为任意值。有 10 个输入字段用于支撑点的 y 值和 x 值。只要阶数小于或等于 10，这些字段就是可编辑的，并且可以操纵支撑点。如果要测试阶数大于 10 的情况，则必须在软件中设置支撑值（我使代码尽可能简单）。支撑点实现在数组 x_k 和 y_k 中，插值波形放入包含“datapoints”元素的数组 xp 和 yp 中。全局常量“deltaX”可能很清楚。但不一定非要使用它。x_k 的值可以自由定义。

波形由函数 DrawGraph(int maxPoints, double minX, double maxX, bool doClear) 绘制在 pGraph 面板中。该函数将图形缩放到适合面板，并在绘制前清除面板（如果 doClear 为 true）。

计算在 button1_Click 事件中完成。

关注点

哪种算法是最好的，确实很大程度上取决于需要插值的波形，因此最好比较这些算法并观察插值形状。对于少量支撑点，简单的多项式算法可能是最佳解决方案。但即使在 10 个支撑点的情况下，也应该使用拉格朗日或牛顿算法，其中牛顿算法更可取。样条算法非常花哨，但并非总是值得使用，因为它实现起来相当复杂，而且调试起来也不太容易（以防万一 :-))。

 

有关更详细的描述，请访问 www.mosismath.com :-)

 

从多项式到自然样条的插值 - CodeProject - 代码之家
© . All rights reserved.