[计算机数值分析]牛顿插值公式

Spring-_-Bear 的 CSDN 博客导航

---

埃特金算法虽然具有承袭性，但其算式是递推型的，不便于进行理论上的分析。所以采用具有承袭性的显式的牛顿插值公式是不错的选择。

$$p_n(x)=f(x_0)+f(x_0,x_1)(x-x_0)+...+f(x_0,x_1,...,x_n,x)(x-x_0)(x-x_1)...(x-x_n)\text{\hspace{1em}(1)}$$

例：已知四个节点

$$(0,1)，(2,3)，(3,2)，(5,5)$$

构造三次牛顿插值多项式，计算当 $x=2.5$ 时牛顿多项式插值结果。

解：根据牛顿插值公式，如式（1），我们很容易写出三次牛顿插值多项式如下：

$$p_3(x)=f(x_0)+f(x_0,x_1)(x-x_0)+f(x_0,x_1,x_2)(x-x_0)(x-x_1)+f(x_0,x_1,x_2,x_3)(x-x_0)(x-x_1)(x-x_2)$$

化简得

$$p_3(x)=\frac{3}{10}{x}^3-\frac{13}{6}{x}^2+\frac{62}{15}x+1$$

当 $x=2.5$ 时，三次牛顿插值多项式结果为 $p_3(x)=2.479167$

运行示例：

程序源码：

#include <iostream>

using namespace std;

#define MAX 10

typedef struct Point
{
    double x;
    double y;
} point;

int main(void)
{
    int n;
    cout << "请输入插值节点的个数（<10）:";
    cin >> n;

    point p[MAX];
    for (int i = 1; i <= n; i++)
    {
        cout << "请输入第 " << i << " 个点的坐标:";
        cin >> p[i].x;
        cin >> p[i].y;
    }

    double x;
    cout << "请输入 x 的值:";
    cin >> x;

    double sum = 0;
    for (int i = 1; i <= n; i++)
    {
        // 根据牛顿插值公式，第一项为 f(x1) 即 0 阶差商
        if (i == 1)
        {
            sum += p[1].y;
            continue;
        }

        // 求右半部分，即 (x-x0) || (x-x0)(x-x1) || (x-x0)(x-x1)(x-x2)...
        double right = 1;
        for (int j = 1; j < i; j++)
        {
            right *= (x - p[j].x);
        }

        // 求左半部分 leftSum 即 f(x0,x1) || f(x0,x1,x2)...也即差商
        double leftSum = 0;
        for (int k = 1; k <= i; k++)
        {
            if (i == 1)
            {
                leftSum = 1;
                break;
            }

            double denominator = 1;
            for (int m = 1; m <= i; m++)
            {
                if (m == k)
                {
                    continue;
                }
                // 累乘分母
                denominator *= (p[k].x - p[m].x);
            }

            // 累加组成差商的每一项，最终结果为差商
            leftSum += p[k].y / denominator;
        }

        // 累加差商与右半部分的和得最终插值结果
        sum += leftSum * right;
    }

    cout << "\n运用牛顿插值法求得 f(x) =" << sum << endl;

    return 0;
}