埃特金算法虽然具有承袭性,但其算式是递推型的,不便于进行理论上的分析。所以采用具有承袭性的显式的牛顿插值公式是不错的选择。
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 ) (1) 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}})\tag1 pn(x)=f(x0)+f(x0,x1)(x−x0)+...+f(x0,x1,...,xn,x)(x−x0)(x−x1)...(x−xn)(1)
例:已知四个节点
( 0 , 1 ) , ( 2 , 3 ) , ( 3 , 2 ) , ( 5 , 5 ) (0,1),(2,3),(3,2),(5,5) (0,1),(2,3),(3,2),(5,5)
构造三次牛顿插值多项式,计算当 x = 2.5 x=2.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)=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}) p3(x)=f(x0)+f(x0,x1)(x−x0)+f(x0,x1,x2)(x−x0)(x−x1)+f(x0,x1,x2,x3)(x−x0)(x−x1)(x−x2)
化简得
p 3 ( x ) = 3 10 x 3 − 13 6 x 2 + 62 15 x + 1 p_{3}(x)=\frac{3}{10}x^{3}-\frac{13}{6}x^{2}+\frac{62}{15}x+1 p3(x)=103x3−613x2+1562x+1
当 x = 2.5 x=2.5 x=2.5 时,三次牛顿插值多项式结果为 p 3 ( x ) = 2.479167 p_{3}(x)=2.479167 p3(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;
}