简介
Lagrange 插值方法虽然易算,但若要增加一个节点时,全部基函数 φ(x) 都需重新算过。本节介绍另外一种插值方法-牛顿插值法。牛顿插值法是一种重要的插值方法,与拉格朗日插值法在同阶时产生的多项式在化简以后是一样的,余项也是一样的,因此两者只是写法形式不同而已。在针对变化的数据时,其各自独特的写法形式给其适合的应用场景带来了区别。牛顿插值法适用于被插的点不断增多,而且可以很好的复用以前的计算结果。
来源
我们知道 Lagrange 插值多项式的插值基函数为:
φ
i
(
x
)
:
=
∏
j
=
0
,
j
≠
i
n
x
−
x
j
x
i
−
x
j
(
i
=
0
,
1
,
2
,
3
,
⋅
⋅
⋅
,
n
)
\ φ_i\left( x \right) :=\prod_{j=0,j\ne i}^n{\frac{x-x_j}{x_i-x_j}} (i=0,1,2,3,···,n)
φi(x):=j=0,j=i∏nxi−xjx−xj(i=0,1,2,3,⋅⋅⋅,n)
上述基函数形式上比较复杂、计算量大,并且重复计算很多,由线性代数可知任何一个 n 次多项式都可以表示成:
1
,
(
x
−
x
0
)
,
(
x
−
x
0
)
(
x
−
x
1
)
,
⋅
⋅
⋅
,
(
x
−
x
0
)
(
x
−
x
1
)
⋅
⋅
⋅
(
x
−
x
n
−
1
)
1 , (x − x_0) , (x − x_0)(x − x_1) , · · · , (x − x_0)(x − x_1) · · · (x − x_{n−1})
1,(x−x0),(x−x0)(x−x1),⋅⋅⋅,(x−x0)(x−x1)⋅⋅⋅(x−xn−1)
共 n + 1 个多项式的线性组合。那么,是否可以将这 n + 1 个多项式作为插值基函数呢?答案是肯定的。多项式组
1
,
(
x
−
x
0
)
,
(
x
−
x
0
)
(
x
−
x
1
)
,
⋅
⋅
⋅
,
(
x
−
x
0
)
(
x
−
x
1
)
⋅
⋅
⋅
(
x
−
x
n
−
1
)
1 , (x − x_0) , (x − x_0)(x − x_1) , · · · , (x − x_0)(x − x_1) · · · (x − x_{n−1})
1,(x−x0),(x−x0)(x−x1),⋅⋅⋅,(x−x0)(x−x1)⋅⋅⋅(x−xn−1)线性无关,因此可以作为插值基函数。
设已知 n + 1 个不同节点 x 0 , x 1 , ⋅ ⋅ ⋅ , x n x_0, x_1, · · · , x_n x0,x1,⋅⋅⋅,xn 上的函数值 f ( x 0 ) = y 0 , f ( x 1 ) = y 1 , ⋅ ⋅ ⋅ , f ( x n ) = y n f(x_0) = y_0, f(x_1) = y_1, · · · , f(x_n) = y_n f(x0)=y0,f(x1)=y1,⋅⋅⋅,f(xn)=yn构造满足插值条件 p n ( x i ) = f ( x i ) = y i ( i = 0 , 1 , ⋅ ⋅ ⋅ , n ) p_n(x_i) = f(x_i) = y_i(i = 0, 1, · · · , n) pn(xi)=f(xi)=yi(i=0,1,⋅⋅⋅,n) 的次数不超过 n 次多项式。
p n ( x ) = a 0 + a 1 ( x − x 0 ) + a 2 ( x − x 0 ) ( x − x 1 ) + a 3 ( x − x 0 ) ( x − x 1 ) ( x − x 2 ) + ⋅ ⋅ ⋅ + a n ( x − x 0 ) ( x − x 1 ) ( x − x 2 ) ⋅ ⋅ ⋅ ( x − x n − 1 ) p_n(x) = a_0 + a_1(x − x_0) + a_2(x − x_0)(x − x_1) + a_3(x − x_0)(x − x_1)(x − x_2) + · · ·+ \\a_n(x − x_0)(x − x_1)(x − x_2)· · ·(x − x_{n−1}) pn(x)=a0+a1(x−x0)+a2(x−x0)(x−x1)+a3(x−x0)(x−x1)(x−x2)+⋅⋅⋅+an(x−x0)(x−x1)(x−x2)⋅⋅⋅(x−xn−1)
为了便于计算系数值
a
0
,
a
1
,
a
2
,
⋅
⋅
⋅
,
a
n
a_0, a_1, a_2, · · · , a_n
a0,a1,a2,⋅⋅⋅,an,我们引入差商的概念。
我们将差商的值作为牛顿插值公式的系数
利用差商的递推定义,可以用递推来计算差商。
如此,我们便得到了牛顿插值公式的基函数以及它的系数,如此我们就可以写出来牛顿插值公式的函数,其 n 次插值多项式为:
p
n
(
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
n
)
(
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_2)(x − x_0)(x − x_1) + · · ·+ \\ f(x_0, x_1, · · · , x_n)(x − x_0)(x − x_1)· · ·(x − x_{n−1})
pn(x)=f(x0)+f(x0,x1)(x−x0)+f(x0,x1,x2)(x−x0)(x−x1)+⋅⋅⋅+f(x0,x1,⋅⋅⋅,xn)(x−x0)(x−x1)⋅⋅⋅(x−xn−1)
实例
例:已知节点 0, 2, 3, 5 对应的函数值为 1, 3, 2, 5,构造三次Newton 插值多项式,当 x = 2.5 时,计算 Newton 多项式插值结果。
解:由前面计算的差商值,我们很容易写出三次 Newton 插值多项式如下:
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
)
=
1
+
1
×
(
x
−
0
)
−
2
3
×
(
x
−
0
)
(
x
−
2
)
+
3
10
×
(
x
−
0
)
(
x
−
2
)
(
x
−
3
)
=
3
10
×
x
3
−
13
6
×
x
2
+
62
15
×
x
+
1
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) \\ = 1 + 1 \times(x − 0) −\ \frac{2}{3}\times(x − 0)(x − 2) + \ \frac{3}{10}\times(x − 0)(x − 2)(x − 3) \\ \\ =\ \frac{3}{10}\times x^3 −\ \frac{13}{6}\times x^2+\ \frac{62}{15}\times x+1
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)=1+1×(x−0)− 32×(x−0)(x−2)+ 103×(x−0)(x−2)(x−3)= 103×x3− 613×x2+ 1562×x+1
当 x = 2.5 时,Newton 多项式插值结果为:
p
3
(
2.5
)
=
3
10
×
2.
5
3
−
13
6
×
2.
5
2
+
62
15
×
2.5
+
1
=
2.479167
p3(2.5) = \ \frac{3}{10}\times2.5^3 −\ \frac{13}{6}\times 2.5^2+\ \frac{62}{15}\times 2.5+1 = 2.479167
p3(2.5)= 103×2.53− 613×2.52+ 1562×2.5+1=2.479167
代码
#include<iostream>
using namespace std;
double n,x[100],y[100],l[100][100] ;
void input()
{
for(int i=0;i<n;i++)
cin>>x[i]>>y[i];
}
double newton(double test)
{
double re = 0;
for(int i=0;i<100;i++) l[i][0]=y[i];
for(int j =1;j<n;j++)
{
for(int i=j;i<n;i++)
{
l[i][j] = (l[i][j-1]-l[i-1][j-1])/(x[i]-x[i-j]);
cout<<i<<" "<<j<<" "<<l[i][j]<<endl;
}
}
double g[100];
for(int i=0;i<100;i++) g[i]=1.0;
for(int i=0;i<n;i++)
{
for(int j=0;j<i;j++)
{
g[i] *= (test -x[j]);
}
re += g[i]*l[i][i];
}
return re;
}
int main()
{
// for(int i=0;i<100;i++) cout<<l[i]<<endl;
cout<<"输入多少个点"<<endl;
cin>>n;
cout<<"输入插值点:"<<endl;
input();
double test;
while(1)
{
cout<<"输入要计算的值或输入 -1 退出"<<endl;
cin>>test;
if(test == -1) break;
double re = newton(test);
cout<<"result: "<<re <<endl;
}
return 0;
}
结果: