Transcript 第五章、插值法和曲线拟合

第五章 插值和曲线拟合
5.1
5.2
5.3
5.4
引言
插值多项式的构造
分段低次插值
最小二乘法
5.1 引言
定义 5.1 设 y= f(x) 在区间[a,b]上连续，在[a,b]内n+1个互不
相同的点 x0 , x1, xn (a  x0    xn  b) 上取值 y0 , y1, yn 。
如果存在一性态较好的简单函数P(x)，使得
P( xi )  yi
(i  0,1,n)
(5 1)
则称P(x)为f(x)的插值函数。这时，我们称[a,b]为插值区间，
称 x0 , x1 , xn 为插值节（结）点，称（5-1）为插值条件，
f(x)为被插函数。求插值函数P(x)的方法称为插值法。
特别地，当P(x)为次数不超过n次的代数多项式时，相
应的插值法称为多项式插值；当P(x)为三角多项式时，相
应的插值法称为三角插值；当P(x)为分段解析函数时，相
应的插值法称为分段插值。其中三角插值主要用于处理周
期函数。本章仅介绍最基本的多项式插值。
定理 5.1 在 n+1 个互异点 x0 , x1 ,, xn 上满足插值条
件 (5-1) 的次数不超过n次的插值多项式 Pn (x) 存在且惟一。
证 记
n
Pn ( x)   ak x k
k 0
即有
( a0 , a1 ,, an 为实数 )
 1 x0
 x0n  a0   y0 

   
n
 x1  a1   y1 
 1 x1
  
      

   
n 
 y 
1 x
a

x
n

n
n

   n
n
(5  2)
(5  2)
i 1
因为 A  Vn ( x0 , x1 , xn )   ( xi  x j )  0
i 1 j 0
所以，解存在且惟一，这说明由式 (5-2) 表示的 Pn (x) 存在
且惟一，证毕。
5.2 插值多项式的构造
5.2.1 拉格朗日插值多项式
对给定的n+1个次数  n次的插值多项式li ( x) (i  0,1,, n)
1 i  j
li ( x j )   ij  
i, j {0,1,n}
0 i  j
(5  4)
(a  x0    xn  b)
解  x0 ,, xi 1, xi 1 n , xn 为li ( x) 的零点
 li ( x)中含有因式 ( x  x j )而此因式已为n次多项式，故应有
j i
j 0
n
li ( x)  Ai  ( x  x j )
j i
j 0
Ai为待定系数
n
再由 1  li ( xi )  Ai  ( xi  x j )
得
j i
j 0
n
 li ( x)  
j i
j 0
(x  x j )
( xi  x j )
Ai 
1
n
 (x  x )
i
j i
j 0
(i  0,1,, n)
j
l0 ( x),l1 ( x),ln ( x) 称为n次拉格朗日（Lagrange）插值基
函数（或称为拉格朗日基本插值多项式）。据之，我们可
构造多项式
n
n
n
i 0
i 0
j i
j 0
Ln ( x)   yi li ( x)   yi 
x  xj
xi  x j
(5  5)
它称为n次拉格朗日插值多项式
n
引进 n+1 多项式函数  n1( x)   ( x  x j )
(5  6)
j 0
和 i ( x)  n1 ( x) /( x  xi )
n
 i ( x)
拉格朗日插值多项式可表示为 Ln ( x)   yi
i ( xi )
i 0
x  x0
x  x1
形为
L
(
x
)

y

y
n=1时称为线性插值，
1
0
1
x0  x1
x1  x0
n=2时称为抛物插值。
定理 5.2 （误差估计定理）
设 f ( n ) ( x) 在[a, b]上连续，f ( n 1) ( x)在(a, b)上存在，
a  x0  x1    xn  b, Ln ( x) 为满足插值条件
Ln ( xi )  yi (i  0,1,  , n)
的 n 次拉格朗日插值多项式，
则存在    ( x)  (min(x0 , x), max(x, xn ))  (a, b)
f ( n 1) ( )
使得 Rn ( x) 
n 1 ( x), x  [a, b]
(n  1)!
(5  7)
这里 Rn ( x)  f ( x)  Ln ( x) 称为用Ln ( x)近似代替f ( x)的
截断误差，或称为Ln ( x)的（插值）余项。
注 （1）余项公式主要用于理论分析。实际使用时，代
之以误差估计式
M n 1
Rn ( x) 
 n 1 ( x)
(n  1)!
（2）插值节点的选取应尽量靠近插值点，以使 n1 ( x)
尽可能小，以减小误差。
推论 若 f ( x) 为次数不超过n 次的多项式，则
f ( x)  Ln ( x)
n
特别地
 l ( x)  1
i 0
i
例 5.1 给定函数表
x
1.2
1.3
1.4
1.5
1.6
1.7
lnx 0.182322 0.262364 0.336472 0.405465 0.470004 0.530628
试分别用线性插值和抛物插值求ln 1.46的近似值并估计误差。
解
作线性插值 得
1.46  1.5
1.46  1.4
ln 1.46 
ln 1.4 
ln 1.5  0.377868
1.4  1.5
1.5  1.4
0.51002
误差为 R1 
(1.46  1.4)(1.46  1.5)  6.12 10 4
2!
作抛物插值 得
(1.46  1.5)(1.46  1.6)
(1.46  1.4)(1.46  1.6)
ln 1.4 
ln 1.5
(1.4  1.5)(1.4  1.6)
(1.5  1.4)(1.5  1.6)
(1.46  1.4)(1.46  1.5)

ln 1.6  0.378402
(1.6  1.4)(1.6  1.5)
0.7289
误差为 R2 
(1.46  1.4)(1.46  1.5)(1.46  1.6)
3!
 4.1105
ln 1.46 
精确值为ln1.46  0.37843643
5.2.2 牛顿均差插值多项式
定义 5.2（均差）
设 f ( x) 在区间[a, b] 上连续, 若 x0  [a, b], 则定义
g1 ( x) 
f ( x)  f ( x0 )
x  x0
( x  x0 ; x  [a, b])
为 f ( x) 在[a, b] 上关于点x0 的一阶均差函数, 记为 f [ x0 , x]; 特别地, 称
f [ x0 , x1 ]  g1 ( x1 ) 
f ( x1 )  f ( x0 )
x1  x0
( x1  x0 ; x1  [a, b])
为 f ( x) 关于点 x0 , x1 的一阶均差。
进一步, 我们定义
g 2 ( x) 
g1 ( x)  g1 ( x1 ) f [ x0 , x]  f [ x0 , x1 ]

x  x1
x  x1
( x  x0 , x1 ; x  [a, b])
为 f ( x) 在 [a, b] 上关于点 x0 , x1 的二阶均差函数, 记为
f [ x0 , x1 , x]; 特别地, 称
f [ x0 , x1 , x2 ]  g 2 ( x2 ) 
f [ x0 , x2 ]  f [ x0 , x1 ]
( x2  x0 , x1 ; x2  [a, b])
x2  x1
为 f ( x) 关于点 x0 , x1, x2 的二阶均差。
一般地，对于[a, b] 上的互异点 x0 ,  , xk 1和 k  1阶
均差函数g k 1 ( x)  f [ x0 ,  , xk  2 , x], 定义
g k 1 ( x)  g k 1 ( xk 1 )
g k ( x) 
x  xk 1

f [ x0 ,  , xk  2 , x]  f [ x0 ,  , xk  2 , xk 1 ]
x  xk 1
( x  x j ; j  0,1,  , k  1; x  [a, b])
为 f ( x) 在 [a, b] 上关于点x0 ,  xk 1 的 k 阶均差函数, 记为
f [ x0 ,  , xk 1 , x]。特别地, 称
f [ x0 ,  , xk 1 , xk ]  g k ( xk )
f [ x0 ,  , xk  2 , xk ]  f [ x0 ,  , xk  2 , xk 1 ]

xk  xk 1
( xk  x j ; j  0,1,  , k  1; xk  [a, b])
为 f ( x) 在点 x0 ,  xk 上的 k 阶均差。
均差又称为差商, 点 x0 , x1 ,  , xk 称为均差的节点。
例 5.2 给定表格函数
0
1
2
3
4
5
6
7
x
f(x) 10 7 -3 -19 -39 -59 -71 -59
试求均差 f [0,4], f [0,2,4], f [0,2,4,8]
解 首先由定义得
8
9
f (4)  f (0)
 12.25
40
f (2)  f (0)  3  10
再由 f [0,2] 

 6.5
20
2
f [0,4]  f [0,2]
得 f [0,2,4] 
 2.875
42
f (8)  f (0)
又由 f [0,8] 
 0.125
80
f [0,8]  f [0,2]
得 f [0,2,8, ] 
 1.0625
82
f [0,2,8]  f [0,2,4]
所以 f [0,2,4,8] 
 0.984375
84
f [0,4] 
均差的性质：
(1)
k
f [ x0 , x1 , , xk ]  
i 0
f ( xi )
(5  12)
k
 (x  x )
i
j
j i
j 0
由此知，均差与节点的排列顺序无关，即有
f [ x0 , x1 , , xk ]  f [ xi0 , xi1 , , xik ]
(i0 , i1 , ik 为 0，
1，
 k 的任一排列）
这性质又称为均差关于自变量对称
（2） 若 Pn ( x) 为n 次多项式，则其k 阶均差函数Pn [ x0 ,, xk 1 , x]
当 k  n时为n  k 次多项式，当k  n时恒为零。
（3） 若 f ( x) 在区间[a, b] 上一阶导数连续,
x j  [a, b]( j  0,1, n)互异，
则定义带重节点的均差
f [ x0 , x0 , x1 , x n ]  lim f [ x, x0 , x1 ,, x n ],则
x  x0
f [ x0 , x0 , x1 ,, x n ]  f ' [ x1 ,, x n , x]
x  x0
例 5.3 试用列表法对例5.2的表格函数求 f[1, 3, 5, 7]
解 列表计算得
xi
f(xi)
1
3
5
7
7
-19
-59
-59
一阶均差 二阶均差
-13
-20
0
所以 f [1, 3, 5, 7] = 1.125
-1.75
5
三阶均差
1.125
由均差定义，我们有
f ( x)  f ( x0 )  f [ x0 , x](x  x0 )
(5  13)
f [ x0 , x]  f [ x0 , x1 ]  f [ x0 , x1 , x](x  x1 )
(5  14)
f [ x0 , x1 , x]  f [ x0 , x1 , x2 ]  f [ x0 , x1 , x2 , x](x  x2 )
(5  15)

f [ x0 ,  , xn 1 , x]  f [ x0 , , xn ]  f [ x0 ,  , xn , x](x  xn )
(5  16)
把（5  14）代入（5  13）得
f ( x)  f ( x0 )  f [ x0 , x1 ](x  x0 )  f [ x0 , x1 , x](x  x0 )(x  x1 )
记为N1 ( x)  R1 ( x)
(5  17)
再把（5  15）代入（5  17）得
f ( x)  f ( x0 )  f [ x0 , x1 ](x  x0 )  f [ x0 , x1 , x2 ](x  x0 )(x  x1 )
 f [ x0 , x1 , x2 , x](x  x0 )(x  x1 )(x  x2 )
 N 2 ( x)  R2 ( x)
这里N 2 ( x)  N1 ( x)  f [ x0 , x1 , x2 ](x  x0 )(x  x1 )
2
2 1
k 1
j 0
 f ( x0 )   f [ x0 , x1 , , xk ] ( x  x j )
一般地，f ( x)  N n ( x)  Rn ( x)
(5  18)
n
N n ( x)  f ( x0 )   f [ x0 , x1 , xk ]k ( x)
k 1
 N n 1 ( x)  f [ x0 , x1 ,, xn ]n ( x)
(5  19)
k 1
其中k ( x)   ( x  x j ) (k  1,2, n)
j 0
Rn ( x)  f [ x0 , x1 ,, xn , x]n 1 ( x)
称为 N n ( x) 近似 f ( x)的余项。
(5  20)
下面证明N n (x) 是一个插值多项式。
证  Ln ( xi )  f ( xi ) (i  0,1,, n)
 Ln [ x0 , x1 ,  , x k ]  f [ x0 , x1 ,  , x k ] (k  1,2,  , n)
又 Ln [ x0 , x1 ,  , x n , x]  0
n
故 Ln ( x)  Ln ( x0 )   Ln [ x0 , x1 ,  , x k ] k ( x)  Ln [ x0 , x1 ,  , x n , x] n 1 ( x)
k 1
n
 f ( x0 )   f [ x0 , x1  , x k ] k ( x)  0
k 1
 N n ( x)
Nn(x)称为牛顿均差插值多项式。
f ( n 1) ( ( x))
进而, 由
n 1 ( x)  f [ x0 , x1 ,, xn , x]n1 ( x)
(n  1)!
f ( n 1) ( ( xn 1 ))
得 f [ x0 , x1 ,, xn , xn 1 ] 
( ( xn 1 )  ( x0 , xn 1 ))
(n  1)!
f ( n ) ( )
从而，f [ x0 , x1 , xn ] 
(  ( x0 , xn ))
(5  21)
n!
实际用N n ( x) 近似 f ( x) 时，我们可选取适当的节点xn 1近似地
估计误差，即有
~
Rn ( x)  f [ x0 , x1 ,  , xn , xn 1 ]n 1 ( x)记为 Rn ( x)
(5  22)
此外，我们可利用插值多项式进行“反插”计算。
例 5.4 给定表格函数
x
1
2
3
4
5
f(x)
0.5
0.175 1.31 -1.495 10.36
（1）试用二次牛顿均差插值法求 f (2.8) 的近似值；
（2）设 f (x)=-1.166 已知，试用（1）中构造的插值多项
式求 x 的近似值。
解 (1) N 2 (2.8)  0.175  1.135(2.8  2)  1.97(2.8  2)(2.8  3)
 1.3982
f (2.8)  N 2 (2.8)  f [2,3,4,1](2.8  2)(2.8  3)(2.8  4)
 0.9  0.8  (0.2)  (1.2)
 0.1728
(2) 由0.175  1.135
（x  2） 1.97
（x  2）
( x  3)  1.166
整理得  1.97x 2  10.985x  12.749  0
解得 x1  3.929，x 2  1.647 （超出区间[2，
4]，舍去）
故所求解为x  3.929
5.3 分段低次插值（略）
5.4 最小二乘法
最小二乘法的提出
数据的多项式最小二乘拟合
5.4.1
5.4.2
给定 n 组数据( xi , yi ) (i  1,2,  , n) , 求 m 次多项式
m
P ( x)   ak* x k
( m  n)
(5  27)
使 F (a0* , a1* ,  , am* )  min F (a0 , a1 ,  , am )
(5  28)
*
m
k 0
n
m
i 1
k 0
这里 F (a0 , a1 ,  , am )   ( ak xik  yi ) 2 称为 方差。
解
由驻点方程
n
m
F
  2[ ak xik  yi ]xil  0 (l  0,1,2,  , m)
al i 1 k 0
m
整理得
n
 x
k  0 i 1
k l
i
n
ak   xil yi (l  0,1,  , m)
i 1
(5  29)

 na0
 n

即 ( xi )a 0
 in1

m
(
x
  i )a 0
 i 1
引进记号


n
( xi )a1
  
( x ) a m

( x )a1


 ( xim 1 )a m

i 1

i 1
n
i 1
n
2
i
 ( xim 1 )a1
i 1
n
  
i 1
n

k
S

x

k
i

i 1

n
Tk   xik yi

i 1
可简写成
 S0

 S1


S
 m
n
S1
S2

S m 1
n
m
i
( xi2 m )a m
i 1

n
y
i 1
n
i
x y
i 1
n
i
m
x
 i yi
i 1
(k  0,1,,2m)
(5  30)
(k  0,1,, m)
S m  a0   T0 
   
 S m 1  a1   T1 
 



 


   
 S 2 m  am   Tm 

它称为法方程组（或正规方程组）。
i
(5  31)
例 5.5 试对以下数据进行多项式拟合
xi
yi
1
2
3
4
5
6
7
1.1 3.8 8.7 15.6 24.6 37.4 49.6
8
64.2
解 令 P( x)  a0  a1 x  a 2 x 2 , 计算得
S 0  8 , S1  36 , S 2  204 , S 3  1296, S 4  8772,
T0  205, T1  1305.4 , T2  8844.8
从而建立法方程组
36
204  a0   205 
 8

  

 36 204 1296 a1   1305.4 
 204 1296 8772 a   8844.8 

 2  

解得 a0  0.1143, a1  0.0548, a2  1.0190
故 P( x)  0.1143 0.0548x  1.0190 x 2
当拟合多项式中x 的某些幂次不出现时，相应的法方程组
即等于从原法方程（5  31）去掉相应的行列后所得到的方程组。
例如，当拟合函数为Q3 ( x)  a1* x  a3* x 3 时，可从不缺项的法方
程组
 S0

 S1
S
 2
S
 3
S1
S2
S3
S4
S2
S3
S4
S5
S3  a0   T0 
   
S 4  a1   T1 
 



S 5 a2
T
   2 
S 6  a3   T3 
去掉第一，三行和一，三列得到方程组
 S2

 S4
S 4  a1   T1 
    
S 6  a3   T3 
(5  32)
当拟合函数为 Q2 ( x)  a0*  a2* x 2 时，相应的法方程组为
 S0

 S1
S
 2
S1
S2
S3
S 2  a0   T0 
   
S3  a1    T1 
S 4  a2   T2 
 S0

 S2
S 2  a0   T0 
    
S 4  a2   T2 
即
(5  33)
~
当拟合函数为 Q2  a* x  b* x 2 时，相应的法方程组为
 S2

 S3
S3  a   T1 
    
S 4  b   T2 
(5  34)
例 5.6 用最小二乘法求形如 y = ax + bx2 的多项式，使与下
列数据拟合（得数保留三位小数）
x
y
-2
1.9
0
0
1
5.1
2
13.8
解  S 2  9，S 3  1 ，S 4  33 ，T1  28.9 ，T2  67.9
 9 1  a   28.9 
   

 相应的法方程组为
 1 33 b   67.9 
解得 a  2.993，b  1.967
故所求拟合多项式为 y  2.993x  1.967x 2
例 5.7 给定数据
x
-5
y
52.6
0
3.4
1
5.5
2
10.5
试求形如 y = a + bx2 的拟合多项式（得数保留三位小数）。
解  S 0  4 ，S 2  30 ，S 4  642，T0  72 ，T2  1362.5，
 相应的法方程组为
 4 30  a   72 

   

 30 642 b  1362.5 
解得 a  3.207，b  1.972
因此，所求拟合多项式为 y  3.207  1.972x 2
5.4.3
最小二乘法的应用例
5.4.3.1 数据的指数拟合
例 5.8 试对下表数据进行指数拟合
xi
1
2
3
4
5
yi
15.4 20.4 26.8 37.2 48.8
6
63.3
解 首先根据 yi 的值算出zi  ln yi（i  1，
2，
，
6），列出下表。
xi
1
2
Zi=ln yi 2.7344 3.0155
3
4
3.2884 3.6163
5
6
3.8877
4.1479
令 z ( x)  a0  a1 x, 对数据（xi , z i ) (i  1,2,,6) 进行线性拟合
 6 21 a0   20.6902

   

建立法方程组
 21 91 a1   77.4217
解出 a0  2.44717, a1  0.286061
所以 z ( x)  2.44717 0.28606x
y ( x)  e z ( x )  11.55556e 0.28606 x
5.4.3.2 超定方程组的最小二乘解
Ax  b 的法方程的形式为
（AT A) x  AT b
(5  39)
例 5.9 试求以下超定方程组的最小二乘解。
 1.2 x 2  1.6
 x1
2.1x  2 x
 2.4

1
2

x2
 2
 3x1 
1.5 x1 
x2
 2
解
因为
16.6 9.9  T
15.64
 , A b  

A A  
 9.9 7.44
10.73
16.6 9.9  x1  15.64
   

所以法方程为
 9.9 7.44 x 2  10.72
解得 x1  0.4014，x 2  0.9067。
T