Transcript 第二章

第2章
插 值 法
第1节 引言
第2节 拉格朗日插值
第3节 均差与牛顿插值多项式
第4节 埃尔米特插值
第5节 分段低次插值
第6节 三次样条插值
1
2.1
引
2.1.1
言
插值问题的提出
设函数 y  f (x) 在区间 [a, b]上有定义，且已知在点
a  x0  x1    xn  b 上的值 y0 , y1 , , yn ，若存在一简
单函数 P (x)，使
P( xi )  yi
(i  0,1,, n),
（1.1）
成立，就称 P(x) 为 f (x) 的插值函数，点 x0 , x1 ,, xn称为插
值节点，包含节点的区间 [a, b] 称为插值区间，求插值函数
P (x)的方法称为插值法.
2
若 P (x) 是次数不超过 n的代数多项式，即
P( x)  a0  a1 x    an x n
（1.2）
其中 ai为实数，就称P ( x )为插值多项式，相应的插值
法称为多项式插值.
若 P ( x )为分段的多项式，就称为分段插值.
若 P ( x )为三角多项式 ，就称为三角插值.
本章只讨论多项式插值与分段插值.
3
从几何上看，插值法就是确定曲线 y  P(x) ，使其通过
给定的 n  1 个点 ( xi , yi ), i  0,1,, n ，并用它近似已知曲线
y  f (x) . 见图.
P(x)  f(x)
x0
4
x1
x2
x
x3
x4
2.1.2
多项式插值
设在区间 [a, b]上给定 n  1 个点
a  x0  x1    xn  b
上的函数值 yi  f ( xi )(i  0,1,, n),求次数不超过 n 的多项式
P(x) ，使
P( xi )  yi
(i  0,1,, n),
（1.3）
由此可以得到关于系数 a0 , a1 , , an 的 n  1 元线性方程组
5
a0  a1 x0    an x0n  y0 ,

n
a0  a1 x1    an x1  y1 ,


a  a x    a x n  y ,
1 n
n n
n
 0
（1.4）
此方程组的系数矩阵为
1

1
A


1
x0
x1



xn

x0n 
n
x1 
,

n
xn 
（1.5）
称为范德蒙德（Vandermonde）矩阵,由于 xi (i  0,1,, n)
互异，故
6
det A 
n 1
 (x
i
 x j )  0.
i , j o
i j
因此线性方程组（1.4）的解 a0 , a1 , , an 存在且唯一.
定理1
唯一的.
7
满足条件（1.3）的插值多项式 P (x )是存在
2.2
2.2.1
拉格朗日插值
线性插值与抛物插值
对给定的插值点，可以用多种不同的方法求得形如(1.2)
的插值多项式.
先讨论 n  1的简单情形.
（1.2）
P( x)  a0  a1 x    an x n
问题： 给定区间 [ xk , xk 1 ] 及端点函数值
yk  f ( xk ), yk 1  f ( xk 1 ) ，
要求线性插值多项式 L1 ( x) ，使它满足
L1 ( xk )  yk ,
L1 ( xk 1 )  yk 1.
其几何意义就是通过两点( xk , yk ), ( xk 1 , yk 1 ) 的直线. 如图2-2.
9
图2-2
由 L1 ( x) 的几何意义可得到表达式
L1 ( x)  yk 
yk 1  yk
( x  xk )
xk 1  xk
xk 1  x
x  xk
L1 ( x) 
yk 
yk 1
xk 1  xk
xk 1  xk
（点斜式），
（2.1）
（两点式），
由两点式看出，L1 ( x) 是由两个线性函数
lk ( x ) 
x  xk 1
,
xk  xk 1
lk 1 ( x) 
x  xk
,
xk 1  xk
（2.2）
的线性组合得到，其系数分别为 y k 及 yk 1 ，即
L1 ( x)  yk lk ( x)  yk 1lk 1 ( x)
1
0
（2.3）
显然，lk (x) 及 lk 1 ( x) 也是线性插值多项式，在节点 xk 及 xk 1
上满足条件
lk ( xk ) 1,
lk ( xk 1 )  0;
lk 1 ( xk )  0,
lk 1 ( xk 1 )  1,
称 lk (x) 及 lk 1 ( x) 为线性插值基函数，图形见图2-3.
1
1
图2-3
1
2
下面讨论 n  2 的情形.
假定插值节点为 xk 1 ，xk ，xk 1 ，要求二次插值多项式
L2 ( x), 使它满足
L2 ( x j )  y j
( j  k 1, k , k 1)
几何上 L2 ( x)是通过三点 ( xk 1 , yk 1 ), ( xk , yk ), ( xk 1 , yk 1 ) 的抛物线.
可以用基函数的方法求 L2 ( x) 的表达式，此时基函数
lk 1 ( x), lk (x), lk 1 ( x) 是二次函数，且在节点上满足条件
lk 1 ( xk 1 )  1,
lk ( xk )  1,
1
3
lk 1 ( x j )  0,
lk ( x j )  0,
lk 1 ( xk 1 )  1,
lk 1 ( x j )  0,
( j  k , k  1);
( j  k  1, k  1);
( j  k  1, k ).
（2.4）
接下来讨论满足（2.4）的插值基函数的求法，
以求 lk 1 ( x)为例，由插值条件，它应有两个零点 xk 及 xk 1，
可表示为
lk 1 ( x)  A( x  xk )( x  xk 1 ),
lk l1 ( x(k x1 ) ) 1, 1, lk l1 ( x(j x) ) 0, 0, ( j(j k, kk,k1);1);
k 1
k 1
k 1
j
lk 1 ( xk 1 )  1 定出
可由插值条件
其中 A为待定系数，
lk ( xk )  1, lk ( x j )  0,
( j  k  1, k  1);
（2.4）
1
A

lk 1 ( xk 1 ) (1x, k 1lk1x(kx)(j )x
0, ( j  k  1, k ).
k 1  xk 1 )
于是
lk 1 ( x) 
1
4
( x  xk )( x  xk 1 )
.
( xk 1  xk )( xk 1  xk 1 )
同理
lk ( x ) 
( x  xk 1 )( x  xk 1 )
.
( xk  xk 1 )( xk  xk 1 )
lk 1 ( x) 
( x  xk 1 )( x  xk )
.
( xk 1  xk 1 )( xk 1  xk )
二次插值基函数 lk 1 ( x) ，lk (x)，lk 1 ( x) 在区间 [ xk 1 , xk 1 ] 上的
图形见图2-4.
1
5
图2-4
1
6
利用 lk 1 ( x)，lk (x)，lk 1 ( x) ，立即得到二次插值多项式
L2 ( x)  yk 1lk 1 ( x)  yk lk ( x)  yk 1lk 1 ( x)
显然，它满足条件 L2 ( x j )  y j , ( j  k  1, k , k  1).
将 lk 1 ( x),lk (x) ，lk 1 ( x) 代入 (2.5) , 得
L2 ( x)  yk 1
( x  xk )( x  xk 1 )
( xk 1  xk )( xk 1  xk 1 )
( x  xk 1 )( x  xk 1 )
 yk
( xk  xk 1 )( xk  xk 1 )
1
7
( x  xk 1 )( x  xk )
 yk 1
.
( xk 1  xk 1 )( xk 1  xk )
（2.5）
2.2.2
拉格朗日插值多项式
将前面的方法推广到一般情形，讨论如何构造通过
n  1 个节点 x0  x1    xn 的 n 次插值多项式 Ln (x) .
根据插值的定义 Ln (x) 应满足
Ln ( x j )  y j
( j  0,1,, n).
为构造 Ln (x)，先定义 n 次插值基函数.
1
8
（2.6）
定义1
若 n次多项式 l j ( x) ( j  0,1,, n)在 n  1个节点
x0  x1    xn 上满足条件
1, k  j;
l j ( xk )  
( j , k  0,1,, n)
0, k  j.
（2.7）
就称这 n  1 个 n次多项式 l0 ( x), l1 ( x), , ln ( x) 为节点
x0 , x1 ,, xn 上的 n 次插值基函数.
1
9
与前面的推导类似，n 次插值基函数为
( x  x0 )  ( x  xk 1 )( x  xk 1 )  ( x  xn )
lk ( x ) 
( xk  x0 )  ( xk  xk 1 )( xk  xk 1 )  ( xk  xn )
(k  0,1,, n).
（2.8）
显然它满足条件（2.7）.
于是，满足条件（2.6）的插值多项式 Ln (x)可表示为
1, k  j;n
l j ( xk )  
( j , k  0,1,, n)
Ln(
0,x )k 
j. y k lk ( x ).
Ln ( x j ) k y0j ( j  0,1,, n).
2
0
（2.7）
（2.9）
（2.6）
由 lk (x) 的定义，知
n
Ln ( x j )   yk lk ( x j )  y j
( j  0,1, , n).
k 0
形如（2.9）的插值多项式Ln (x)称为拉格朗日插值多项式，
而(2.3)与(2.5)是 nn  1 和 n  2 的特殊情形.
（2.9）
Ln ( x )   y k lk ( x ).
若引入记号 k  0
n1 ( x)  ( x  x0 )( x  x1 )( x  xn ),
L1 ( x)  yk lk ( x)  yk 1lk 1 ( x)
容易求得
（2.10）
（2.3）
（2.5）

n1 ( xk )  ( xk  x0 )( xk  xk 1 )( xk  xk 1 )( xk  xn ),
L2 ( x)  yk 1lk 1 ( x)  yk lk ( x)  yk 1lk 1 ( x)
2
1
于是公式（2.9）可改写成
n 1 ( x)
Ln ( x)   yk
.

(
x

x
)

(
x
)
k 0
k
n 1
k
n
（2.11）
注意: n 次插值多项式 Ln (x) 通常是次数为 n 的多项式，
n
Ln ( x )  
n .y k lk ( x ).
特殊情况下次数可能小于
k 0
（2.9）
例如通过三点 ( x0 , y0 ), ( x1 , y1 ), ( x2 , y2 ) 的二次插值多项式
L2 ( x)，如果三点共线，则 y  L2 ( x) 就是一条直线，而不是
抛物线，这时 L2 ( x)是一次多项式.
2
2
n1 ( x)  ( x  x0 )( x  x1 )( x  xn ),
2.2.3
插值余项与误差估计
（2.10）
若在 [ a, b] 上用 Ln (x) 近似 f (x) ， 则其截断误差为
Rn ( x)  f ( x)  Ln ( x), 也称为插值多项式的余项.
定理2
设 f ( n ) ( x) 在[a, b]上连续， f ( n1) ( x)在 (a, b) 内
存在，节点 a  x0  x1    xn  b, Ln ( x)是满足条件(2.6)
的插值多项式，则对任何 x  [a, b] ，插值余项
f ( n 1 ( )
Rn ( x)  f ( x)  Ln ( x) 
n 1 ( x)
Ln ( x j )  y j ( j  0(n,1,
,
1)! n).
（2.14）
（2.6）
这里   (a, b)且依赖于 x ，n1 ( x) 是（2.10）所定义的.
2
3
证明
由给定条件知 Rn (x) 在节点 xk (k  0,1,, n)
上为零，即 Rn ( xk )  0 (k  0,1,, n)，于是
Rn ( x)  K ( x)( x  x0 )( x  x1 )( x  xn )  K ( x)n1 ( x)
（2.13）
其中 K (x) 是与 x有关的待定函数.
现把 x 看成 [a, b]上的一个固定点，作函数
 (t )  f (t )  Ln (t )  K ( x)(t  x0 )(t  x1 ) (t  xn ),
根据 f 的假设可知  ( n ) (t )在 [ a, b]上连续， ( n1) (t ) 在 (a, b)
内存在.
2
4
根据插值条件及余项定义，可知  (t ) 在点 x0 , x1 ,, xn 及
x 处均为零，故  (t ) 在 [a, b] 上有 n  2个零点，
根据罗尔定理，  (t )在  (t )的两个零点间至少有一个零点，
故  (t ) 在 [a, b]内至少有 n  1个零点.
对  (t ) 再应用罗尔定理，可知  (t ) 在 [a, b]内至少有
n 个零点.
依此类推， ( n1) (t )在 (a, b)内至少有一个零点，记为
  (a, b) ，使
 ( n1) ( )  f ( n1) ( )  (n  1)! K ( x)  0,
2
5
于是
f ( n 1) ( )
K ( x) 
,   (a, b),
(n  1)!
且依赖于 x
将它代入（2.13），就得到余项表达式（2.12）.
余项表达式只有在 f (x) 的高阶导数存在时才能应用.
Rn ( x()a, bK)内的具体位置通常不可能给出，
( x)( x  x0 )( x  x1 )( x  xn )  K ( x)n1 ( x) （2.13）
但在
如果可以求出 max f ( n 1) ( x)  M n 1 , 那么插值多项式 Ln (x)
a  x b
逼近 f (x) 的截断误差限是
Rn ( x) 
2
6
M n 1
n 1 ( x) .
(n  1)!
（2.14）
当 n  1 时，线性插值余项为
R1 ( x) 
1
1
f ( )2 ( x)  f ( )( x  x0 )( x  x1 ),
2
2
  [ x0 , x1 ]
（2.15）
当 n  2 时，抛物插值余项为
1
R2 ( x)  f ( )( x  x0 )( x  x1 )( x  x2 ),
6
 [ x0 , x2 ]
2
7
（2.16）
利用余项表达式（2.12），当 f ( x)  x k (k  n)时，
由于 f ( n 1) ( x)  0 ,于是有
n
Rn ( x)  x k   xik li ( x)  0,
i 0
由此得
n
k
k
x
l
(
x
)

x
, k  0,1, , n.
 ii
（2.17）
i 0
特别当 k  0 时，有
n
 l ( x)  1.
i 0
2
8
i
（2.18）
利用余项表达式（2.12）还可知，若被插函数 f ( x)  H n
由于 f ( n 1) ( x)  0 ，故 Rn ( x)  f ( x)  Ln ( x)  0 ，即它的插
值多项式 Ln ( x)  f ( x).
2
9
例1
5
证明  ( xi  x) 2 li ( x)  0 ，其中 li (x) 是关于点
i 0
x0 , x1 ,, x5 的插值基函数.
证明
利用公式（2.17）可得
5
5
2
2
(
x

x
)
l
(
x
)

(
x

2
x
x

x
)li ( x)
 i
 i
i
i
2
i 0
i 0
5
5
5
i 0
i 0
i 0
  xi 2li ( x)  2x xi li ( x)  x 2 li ( x)
 x 2  2 x 2  x 2  0.
3
0
例2 已知 sin 0.32  0.314567, sin 0.34  0.333487,
sin 0.36  0.352274, 用线性插值及抛物插值计算 sin 0.3367
的值并估计截断误差.
解
由题意, 取
y
 yk ,
xL0 1( x0).32
 ,yk y0  k0.1314567
( x  xk )
xk 1  xk
x1  0.34,
y1  0.333487,
x2  0.36,
y2  0.352274.
（点斜式），
用线性插值计算，取 x0  0.32, x1  0.34, 由公式（2.1）
3
1
sin 0.3367  L1 (0.3367)
y1  y0
 y
(0.3367  x0 )
x1  x0
 0.314567 
 0.330365.
3
2
0.01892
 0.0167
0.02
由（2.15）,其截断误差
M2
( x  x0 )( x  x1 ) ,
2
1
1


R1 ( x)  f ( )2 ( x)  f ( )( x  x0 )( x  x1 ),
R1 ( x) 
其中
2
2
M 2  max f ( x)  max  sin x  sin x1  0.3335,
x0  x  x1
x0  x  x1
于是
R1 (0.3367)  sin 0.3367  L1 (0.3367)
1
  0.3335  0.0167  0.0033
2
3
3
 0.92 105.
用抛物插值计算，由公式（2.5）得
sin 0.3367  y0
( x  x0 )( x  x2 )
( x  x1 )( x  x2 )
 y1
( x0  x1 )( x0  x2 )
( x1  x0 )( x1  x2 )
( x  x0 )( x  x1 )
y
L2 ( x)  yk 1lk 1 (2x)( x2 ykxl0k)(( x2) xy1 )k 1lk 1 ( x) （2.5）
 L2 (0.3367)
0.7689 10 4
 0.314567 
 0.333487
0.0008
3.89 10 4
 0.551110 4

 0.352274 
0.0004
0.0008
 0.330374
3
4
这个结果与6位有效数字的正弦函数表完全一样，
这说明查表时用二次插值精度已相当高了.
由（2.14）, 截断误差限
M3
R2 ( x) 
( x  x0 )( x  x1 )( x  x2 ) ,
6
其中
M 3  max f ( x)  cos x0  0.9493,
x0  x  x2
于是
3
5
R2 (0.3367)  sin 0.3367  L2 (0.3367)
1
  0.9493  0.0167  0.033  0.0233
6
 2.0132 106.
3
6
例2
2
设 f  C [a, b] ，试证
max f ( x)  [ f (a) 
a  x b
f (b)  f (a)
1
( x  a)  (b  a) 2 M 2 ,
ba
8
其中 M 2  max f ( x) .
a  x b
证明
通过两点 (a, f (a)) 及 (b, f (b)) 的线性插值为
L1 ( x)  f (a) 
于是
max f ( x)  [ f (a) 
a  x b
3
7
f (b)  f (a)
( x  a),
ba
f (b)  f (a)
( x  a)]
ba
f ( )
 max f ( x)  L1 ( x)  max
( x  a)( x  b)
a  x b
a  x b
2
M
1
 2 max ( x  a)( x  b)  (b  a)2 M 2 .
2 a  x b
8
3
8
3
9
Matlab 实现
 例 对正弦曲线上的数据点(0, 0),
(3pi/2, -1)进行多项式插值
4
0
(pi/2, 1),
(pi, 0),
4
1
4
2
4
3
 Polyval还可以计算多项式在向量上的值
4
4
4
5
4
6
 也可以使用内部命令polyfit进行插值
4
7
 一维插值函数interp1()
 y1=interp1(x,y,x1,方法)
 其中x,y是插值数据，x1为用户指定的一组新的插值点的




横坐标，可以是标量、向量或矩阵。方法：
默认为‘linear’（线性插值，在两个样本点间简单的采
用直线拟合）
‘nearest’（最近点等值方式）
‘cubic’（三次Hermite插值，新版本改为‘pchip’）
‘spline’（三次分段样条插值，建议使用,端点处的信息系
统自动选取）
4
8
 【例】已知的数据点来自函数
 根据生成的数据进行插值处理，得出较平滑的曲线。
 根据给出的函数可以直接生成数据，并绘图
4
9
5
0
5
1
5
2
5
3
5
4
2.3
2.3.1
均差与牛顿插值公式
插值多项式的逐次生成
利用插值基函数很容易得到拉格朗日插值多项式，公
式结构紧凑，在理论分析中甚为方便，但当插值节点增减
时全部插值基函数 lk ( x)( k  0,1,, n)均要随之变化，整个
公式也将发生变化，甚为不便.为了计算方便可重新设计一
种逐次生成插值多项式的方法.
当 n  1时，记线性插值多项式为 P1 ( x) ，插值条件为
P1 ( x0 )  f ( x0 ), P1 ( x1 )  f ( x1 ), 由点斜式
P1 ( x)  f (a) 
f (b)  f (a )
( x  x0 ),
ba
将 P1 ( x)看成是零次插值 P1 ( x)  f ( x0 )的修正，即
P1 ( x)  P0 ( x)  a1 ( x  x0 ),
其中 a1  f ( x1 )  f ( x0 ) 是函数 f (x) 的差商.
x1  x0
对于三个节点的二次插值 P2 ( x) ，插值条件为
P2 ( x0 )  f ( x0 ), P2 ( x1 )  f ( x1 ), P2 ( x2 )  f ( x2 ),
5
6
插值多项式
P2 ( x)  P0 ( x)  a2 ( x  x0 )( x  x1 ).
显然
P2 ( x0 )  f ( x0 ), P2 ( x1 )  f ( x1 ),
由 P2 ( x2 )  f ( x2 ), 得
f ( x2 )  f ( x0 ) f ( x1 )  f ( x0 )

x2  x0
x1  x0
P ( x )  P1 ( x2 )
a2  2 2

.
( x2  x0 )( x2  x1 )
x2  x1
系数 a2是函数 f 的“差商的差商”.
5
7
一般情况，已知 f 在插值点 xi (i  0,1,, n) 上的值
为 f ( xi )(i  0,1,, n)，要求 n次插值多项式 Pn (x)满足
条件
Pn ( xi )  f ( xi ),
i  0,1,, n,
（3.1）
则 Pn (x)可表示为
Pn ( x)  a0  a1 ( x  x0 )    an ( x  x0 ) ( x  xn1 ), （3.2）
其中 a0 , a1 ,, an 为待定系数，可由插值条件确定.
这里的 Pn (x)是由基函数 {1, ( x  x0 ), , ( x  x0 )( x  xn1 )}
逐次递推得到的，这一点与拉格朗日插值不同.
5
8
2.3.2
定义2 称
均差及其性质
f [ x0 , xk ] 
f ( xk )  f ( x0 )
xk  x0
于点 x0 , xk 的一阶均差.
f [ x0 , x1 , xk ] 
称为 f (x) 的二阶均差.
5
9
f [ x0 , xk ]  f [ x0 , x1 ]
xk  x1
为函数 f (x)关
一般地，称
f [ x0 , , xk  2 , xk ]  f [ x0 , x1 , xk 1 ]
f [ x0 , x1 ,  xk ] 
xk  xk 1
（3.3）
为 f (x)的 k 阶均差（均差也称为差商）.
6
0
均差有如下的基本性质：
（1） k 阶均差可表为函数值 f ( x0 ), f ( x1 ), , f ( xk ) 的线
性组合，即
f [ x0 , x1 , , xk ] 
f (x j )
k
 (x
j 0
j
 x0 )  ( x j  x j 1 )( x j  x j 1 )  ( x j  xk )
.
（3.4）
这个性质可用归纳法证明.
这性质也表明均差与节点的排列次序无关，称为均差
的对称性.
6
1
即
f [ x0 , x1 , xk ]  f [ x1 , x0 , x2 ,, xk ]    f [ x1 ,, xk , x0 ]
（2） 由性质（1）及（3.3）可得
f [ x1 , , xk ]  f [ x0 , , xk 1 ]
f [ x0 , x1 ,  xk ] 
.
xk  x0
（3.3）’
f [ x , , xk  2n
, x阶导数，且节点
（3）
k ]  f [ x0 , x1  , xk 1 ]
f (x
f [ x , x若
,
x ) ]在
[ a, b]0 上存在
0
1
xk  xk 1
k
x0 , x1 ,, xn [a, b], 则 n阶均差与导数关系如下：（3.3）
f [ x0 , x1 ,  xk ] 
6
2
f
( )
,
n!
(n)
这公式可直接用罗尔定理证明.
  [a, b].
（3.5）
证明：设 q( x )  f [ x，x0，，xn ]( x  x0 )( x  x1 )( x  xn ), q( x )在
x0，，xn处均为零，所以q( x )在[a , b]上有n  1个零点，根据罗尔定理，
q( x )在q( x )的两个零点间至少有一个零点，故q( x )在[a , b]内至少有n个
零点；反复应用罗尔定理，可知q ( n ) ( x )在[a , b]内至少有1 个零点, 记为
  [a , b]，使
q( n) ( )  f ( n) ( )  n! f [ x0，
，xn ]  0，
所以
63
f ( n ) ( )
f [ x0，，xn ] 
.
n!
均差计算可列均差表如下（表2-1）.
表2  1
6
4
xk
x0
x1
x2
f ( xk ) 一阶均差
f ( x0 )
f ( x1 )
f [ x0 , x1 ]
f ( x2 )
f [ x1 , x2 ]
f [ x0 , x1 , x2 ]
x3
x4

f ( x3 )
f ( x4 )

f [ x1 , x2 , x3 ]
f [ x2 , x3 , x4 ]

f [ x2 , x3 ]
f [ x3 , x4 ]

二阶均差
三阶均差
四阶均差
f [ x0 , x1 , x2 , x3 ]
f [ x1 , x2 , x3 , x4 ]

f [ x0 , x1 , x2 , x3 , x4 ]

2.3.3
牛顿插值公式
根据均差定义，一次插值多项式为
P1 ( x)  P0 ( x)  f [ x0 , x1 ]( x  x0 )  f ( x0 )  f [ x0 , x1 ]( x  x0 ),
二次插值多项式为
P2 ( x)  P1 ( x)  f [ x0 , x1 , x2 ]( x  x0 )( x  x1 )
 f ( x0 )  f [ x0 , x1 ]( x  x0 )  f [ x0 , x1 , x2 ]( x  x0 )( x  x1 ).
6
5
根据均差定义，把 x 看成 [a, b]上一点，可得
f ( x)  f ( x0 )  f [ x, x0 ]( x  x0 ),
f [ x, x0 ]  f [ x0 , x1 ]  f [ x, x0 , x1 ]( x  x1 ),

f [ x, x0 ,, xn1 ]  f [ x0 , x1 ,, xn ]  f [ x, x0 , x1 ,, xn ]( x  xn ).
6
6
只要把后一式依次代入前一式，就得到
f ( x)  f ( x0 )  f [ x0 , x1 ]( x  x0 )
 f [ x0 , x1 , x2 ]( x  x0 )( x  x1 )  
 f [ x0 , x1 ,, xn ]( x  x0 ) ( x  xn1 )
 f [ x, x0 ,, xn ]n1 ( x)
 Pn ( x)  Rn ( x),
其中
Pn ( x)  f ( x0 )  f [ x0 , x1 ]( x  x0 )
（3.6）
 f [ x0 , x1 , x2 ]( x  x0 )( x  x1 )  
 f [ x0 , x1 ,, xn ]( x  x0 ) ( x  xn1 ),
6
7
Rn ( x)  f ( x)  Pn ( x)  f [ x, x0 ,, xn ]n1 ( x),
（3.7）
n1 ( x)是由（2.10）定义的.
显然，由（3.6）确定的多项式 Pn (x) 满足插值条件，
且次数不超过 n，它就是形如（3.1）的多项式，其系数为
1,, n).
（2.10）
n1 ( xa)k ( xf [x0x,0
)( x, xk ],
x1 )((xkx0n,),
称 Pn (x) 为牛顿（Newton）均差插值多项式.
Pn ( x)  f ( x0 )  f [ x0 , x1 ]( x  x0 )
（3.6）
Pn ( x)  a0  a1 ( x  x0 )  a2 ( x  x0 )( x  x1 )  
系数 a k 就是均差表2-1中加横线的各阶均差，它比
 f [ x0a, nx(1 ,xx2 ](
 x(0x)(xx
x )   （3.1）
x0x) 
n 1 )1
拉格朗日插值计算量省，且便于程序设计.
 f [ x0 , x1 ,, xn ]( x  x0 ) ( x  xn1 ),
6
8
（3.7）为插值余项，由插值多项式唯一性知，它与
拉格朗日插值多项式的余项应该是等价的.
事实上，利用均差与导数关系式就可以证明这一点.
Rn ( x)但（3.7）更有一般性，它在
 f ( x)  Pn ( x)  f [ x, x0 ,, xn f]是由离散点给出的
（3.7）
n 1 ( x),
情形或 f 导数不存在时也是适用的.
牛顿插值多项式的优点还在于它的递进性，当增加
插值节点时，只要在原来插值多项式的基础上增加一项
即可.
6
9
例4
给出 f (x)的函数表（见表2-2），求4次牛顿插
值多项式，并由此计算 f (0.596) 的近似值.
首先根据给定函数表造出均差表.
表22
xk
f(xk ) 一阶均差 二阶均差 三阶均差 四阶均差 五阶均差
0.40 0.41075
0.55 0.57815 1.11600
0.65 0.69675 1.18600 0.28000
0.80 0.88811 1.27573 0.35893 0.19733
0.90 1.02652 1.38410
1.05 1.25382 1.51533
7
0
0.43348
0.52493
0.21300
0.22863
0.03134
0.03126
0.00012
从均差表看到4阶均差近似常数，5阶均差近似为0.
故取4次插值多项式 P4 ( x) 做近似即可.
按牛顿插值公式，将数据代入
P4 ( x)  0.41075  1.116( x  0.4)  0.28( x  0.4)( x  0.55)
 0.19733( x  0.4)( x  0.55)( x  0.65)
 0.03134( x  0.4)( x  0.55)( x  0.65)( x  0.8),
于是
f (0.596)  P4 (0.596)  0.63192,
7
1
截断误差
R4 ( x)  f [ x0 ,, x5 ]5 (0.596)  3.63 109.
这说明截断误差很小，可忽略不计.
7
2
2.3.4
差分形式的牛顿插值公式
实际应用时经常遇到等距节点，即 xk  x0  kh (k  0,1,, n)
的情形，这里 h为常数，称为步长， 这时插值公式可以进一步
简化，计算也简单得多.
设 xk 点的函数值为 f k  f ( xk )(k  0,1,, n) ，称 f k  f k 1  f k ,
为 xk 处以 h 为步长的一阶（向前）差分.
类似地称 2 f k  f k 1  f k 为 xk处的二阶差分.
一般地称
n f k  n1 f k 1  n1 f k
7
3
为 xk处的 n 阶差分.
（3.8）
为了表示方便，再引入两个常用算子符号
I fk  fk ,
E f k  f k 1 ,
I 称为不变算子 ，E称为步长为 h 的移位算子，由此
 f k  f k 1  f k  E f k  I f k  (E  I) f k ,
 f k  ( E  I) f k
n
n
n
 n  n j
j n
  (1)  E f k   (1)   f n  k  j , （3.9）
j 0
j 0
 j
 j
n
j
 n  n(n  1)  (n  j  1)
其中   
为二项式展开系数，（3.9）
j!
 j
说明各阶差分均可由函数值给出.
7
4
反之，由
n j
 (I  ) f k  [   ] f k ,
j 0  j 
n
f nk  E f k
n
n
可得
n j
    f k .
j 0  j 
n
f nk
（3.10）
从而有均差与差分的关系：
f k 1  f k
f k
f [ xk , xk 1 ] 

,
xk 1  xk
h
1 2
f [ xk 1 , xk  2 ]  f [ xk , xk 1 ]

 fk ,
f [ xk , xk 1 , xk  2 ] 
2
2h
xk  2  x k
7
5
一般地有
f [ x k ,  , xk  m ] 
1 1 m
 fk ,
m
m! h
(m  1,2, , n).
（3.11）
由（3.11）和（3.5）又可得到差分与导数的关系：
n f k  h n f ( n ) ( ),
其中   ( xk , xk  n ) .
7
6
（3.12）
由给定函数表计算差分可由以下形式差分表给出.
fk

2
3
4

f0
f 0
2 f 0
f1
f1
3 f 0
2 f1
f2
f 2
3 f1
2 f 2
f3
f 3

7
7




f4

4 f 0
在牛顿插值公式（3.6）中，用（3.11）的差分代替均差，
并令 x  x0  th ，则得
t (t  1) 2
Pn ( x0  th)  f 0  tf 0 
 f0  
2!
t (t  1)  (t  n  1) n

 f0 ,
n!
（3.13）
（3.13）称为牛顿前插公式，由（3.7）式得其余项为
t (t  1)  (t  n) n 1 ( n 1)
Rn ( x) 
h f
( ),
(n  1)!
7
8
  ( x0 , xn ).
（3.14）
例5
给出 f ( x)  cos x 在 xk  kh, k  0,1,,5, h  0.1
处的函数值，试用4次牛顿前插公式计算 f (0.048) 的近似
值并估计误差.
解
表 23
为使用牛顿插值公式，先构造差分表.
差分表
xk
f ( xk )
0.00
1.00000
f
2 f
3 f
4 f
5 f
 0.00500
0.10
 0.00993
0.99500
 0.01493
0.20
0.00013
 0.00980
0.98007
 0.02473
0.30
 0.00955
0.95534
 0.4348
7
9
0.50
0.87758
0.00010
0.00035
 0.00920
0.92106
 0.00002
0.00025
 0.03428
0.40
0.00012
取 x  0.048, h  0.1,
则 t
x  x0 0.048  0

 0.48, 得
h
0.1
f (0.048)  cos 0.048  P4 (0.048)
 1.00000  0.48  (0.00500)
(0.48)(0.48  1)
(0.00993)
2
1
 (0.48)(0.48  1)(0.48  2)(0.00013)
3!
1
 (0.48)(0.48  1)(0.48  2)(0.48  3)(0.00012)
4!

 0.99885
8
0
由（3.14）可得误差估计
M5
R4 (0.048) 
t (t  1)(t  2)(t  3)(t  4) h 5
5!
 1.5845 107 ,
t (t  1)  (t  n) n 1 ( n 1)
Rn ( x) 其中 M 5  sin 0h.6 f 0.565
( ),
.
(n  1)!
8
1
  ( x0 , xn ).
（3.14）
2.4 埃尔米特插值
有些实际的插值问题不但要求在节点上函数值相等，
而且还要求对应的导数值也相等，甚至要求高阶导数也相等.
满足这种要求的插值多项式就是埃尔米特插值多项式.
2.4.1
重节点均差与泰勒插值
关于均差，有
定理3
设 f  C n [a, b], x0 , x1 ,, xn 为 [a, b] 上的相异节
点，则 f [ x0 , x1 ,, xn ]是其变量的连续函数.
如果 [a, b] 上的节点互异，根据差商定义，若 f  C1[a, b],
则有
lim f [ x0 , x]  lim
x  x0
x  x0
f ( x)  f ( x0 )
 f ( x0 ).
x  x0
由此定义重节点均差
f [ x0 , x0 ]  lim f [ x0 , x]  f ( x0 ).
8
3
x  x0
类似地可定义重节点的二阶均差，当 x1  x0时，有
f [ x0 , x0 , x1 ] 
f [ x0 , x1 ]  f [ x0 , x0 ]
.
x1  x0
当 x1  x0 时，有
f [ x0 , x0 , x0 ]  lim f [ x0 , x1 , x2 ] 
x1  x0
x 2  x0
1
f ( x0 ).
2
一般地，可定义 n 阶重节点的均差，由（3.5）得
f [ x0 , x0 , , x0 ]  lim f [ x0 , x1 , , xn ] 
xi  x0
8
4
1 (n)
f ( x0 ). （4.1）
n!
在牛顿均差插值多项式中若令 xi  x0 (i  1,2,, n) 则
由（4.1）可得泰勒多项式
f ( n ) ( x0 )
Pn ( x)  f ( x0 )  f ( x0 )( x  x0 )   
( x  x0 ) n . （4.2）
n!
这实际上是在点 x0 附近逼近 f (x)的一个带导数的插值多项
式，它满足条件
Pn( k ) ( x0 )  f ( k ) ( x0 ),
k  0,1,, n.
（4.3）
（4.2）称为泰勒插值多项式，它就是一个埃尔米特插值多
项式，余项为
8
5
f ( n 1) ( )
Rn ( x) 
( x  x0 ) n 1 ,
(n  1)!
  (a, b),
（4.4）
它与插值余项（2.12）中令 xi  x0 (i  1,2,, n) 的结果是
一致的.
泰勒插值是牛顿插值的极限形式，是只在一点 x0 给出
n  1 个插值条件所得到的 n 次埃尔米特插值多项式.
一般地只要给出 m  1个插值条件（包括函数和导数值）
就可以构造出次数不超过 m次的埃尔米特多项式.
8
6
2.4.2
两个典型的埃尔米特插值
考虑满足条件 P( xi )  f ( xi ) (i  0,1,2) 及 P( x1 )  f ( x1 )
的插值多项式及其余项表达式.
由给定的4个条件，可确定次数不超过3的插值多项式.
由于此多项式通过点 ( x0 , f ( x0 )), ( x1 , f ( x1 )), ( x2 , f ( x2 )),
故其形式为
P( x)  f ( x0 )  f [ x0 , x1 ]( x  x0 )
 f [ x0 , x1 , x2 ]( x  x0 )( x  x1 )
 A( x  x0 )( x  x1 )( x  x2 ),
8
7
待定常数 A，可由条件 P( x1 )  f ( x1 ) 确定，
通过计算可得
A
f ( x1 )  f [ x0 , x1 ]  ( x1  x0 ) f [ x0 , x1 , x2 ]
.
( x1  x0 )( x1  x2 )
为了求出余项 R( x)  f ( x)  P( x) 的表达式，可设
R( x)  k ( x)( x  x0 )( x  x1 ) 2 ( x  x2 ),
其中 k (x) 为待定函数.
8
8
构造
 (t )  f (t )  P(t )  k ( x)(t  x0 )(t  x1 ) 2 (t  x2 ).
显然  ( x j )  0 ( j  0,1,2), 且  ( x1 )  0,  ( x)  0,
故  (t )在 (a, b) 内有5个零点（二重根算两个）.
反复应用罗尔定理，得  ( 4) (t ) 在 (a, b) 内至少有一个
零点ξ， 故有
 ( 4) ( )  f ( 4) ( )  4!k ( x)  0,
8
9
于是
k ( x) 
1
f
4!
( 4)
( ),
余项表达式为
1 ( 4)
R( x)  f ( )( x  x0 )( x  x1 ) 2 ( x  x2 ),
4!
式中  位于 x0 , x1 , x2 和 x 所界定的范围内.
9
0
(4.5)
例6
给定 f ( x)  x 3 / 2 , x0  1 , x1  1, x2  9 , 试求 f (x)
4
4
1 9
在 [ , ]上的三次埃尔米特插值多项式 P (x ) ，使它满足
4 4
P( xi )  f ( xi ) (i  0,1,2), P( x1 )  f ( x1 ),
并写出余项表达式.
解
由所给节点可求出
1
1
9
27
f 0  f ( )  , f1  f (1)  1, f 2  f ( ) 
,
4
8
4
8
f ( x) 
3 1/ 2
x ,
2
f (1) 
为了构造牛顿均差插值，先构造均差表
9
1
3
.
2
于是有
xi
fi
1
4
1
8
1
1
9
1 7
1
11
1
P( x)   ( x  ) 
( x  )( x  1) 4
8 6
4
30
4
1
9
 A( x  )( x  1)( x  ).
4
4
再由条件 P(1)  f (1)  3 ，可得
2
7 11 3
3
5
3
P(1)  
 A ( )  ,
6 30 4
4
4
2
27
8
7
11
f [ x0 , x1 ]  , f [ x0 , x1 , x2 ] 
.
6
30
令
9
2
表2 - 4
均差表
7
6
19
10
11
30
解出
A
16 3 7 11
14
(  
)
.
15 2 6 40
225
于是所求的三次埃尔米特多项式为
1 7
1
11
1
 (x  ) 
( x  )( x  1)
8 6
4
30
4
14
1
9

( x  )( x  1)( x  )
225
4
4
14 3 263 2 233
1

x 
x 
x
,
225
450
450
25
P( x) 
9
3
余项
f ( 4 ) ( )
1
9
R( x)  f ( x)  P( x) 
( x  )( x  1) 2 ( x  )
4!
4
4
1 9 5 / 2
1
9
1 9


( x  )( x  1) 2 ( x  )   ( , ).
4! 16
4
4
4 4
9
4
另一个典型例子是两点三次埃尔米特插值，插值节点
取为 xk及 xk 1，插值多项式为 H 3 ( x)，插值条件为
H 3 ( xk )  yk ,
H 3 ( xk )  mk ,
H 3 ( xk 1 )  yk 1 ;



H 3 ( xk 1 )  mk 1. 

(4.6)
采用基函数的方法，令
H 3 ( x)   k ( x) yk   k 1 ( x) yk 1   k ( x)mk   k 1 ( x)mk 1 ,
(4.7)
其中  k ( x),  k 1 ( x),  k ( x),  k 1 ( x) 是关于节点 xk 及 xk 1的
三次埃尔米特插值基函数，分别满足
9
5
 k ( xk )  1,  k ( xk 1 )  0,
 k ( xk )   k ( xk 1 )  0,
 k 1 ( xk )  0,
 k 1 ( xk 1 )  1,
 k 1 ( xk )   k 1 ( xk 1 )  0;
 k ( xk )   k ( xk 1 )  0,
 k ( xk )  1,
 k ( xk 1 )  0,
 k 1 ( xk )   k 1 ( xk 1 )  0,
 k1 ( xk )  0,
9
6
 k1 ( xk 1 )  1.
根据给定条件可令
2
 x  xk 1 
 ,
 k ( x)  ax  b 
 xk  xk 1 
显然
 k ( xk 1 )  0,  k ( xk 1 )  0,
再利用
 k ( xk )  axk  b  1,
及
 k ( xk )  2
9
7
axk  b
 a  0,
xk  xk 1
解得
a
2
,
xk  xk 1
b  1
2 xk
,
xk  xk 1
于是求得

x  xk

 k ( x )  1  2
xk 1  xk

2
 x  xk 1 

 .
 xk  xk 1 
(4.8)
同理可得

x  xk 1  x  xk


 k 1 ( x)  1  2
xk  xk 1  xk 1  xk

9
8



2
(4.9)
为求  k (x) ，由给定条件可令
2
 x  xk 1 
 ,
 k ( x)  a x  xk 
 xk  xk 1 
直接由  k ( xk )  a  1 ，得到
2
 x  xk 1 
 ,
 k ( x)   x  xk 
 xk  xk 1 
(4.10)
同理
 x  xk
 k 1 ( x)   x  xk 1 
 xk 1  xk
9
9
2

 .

(4.11)
最后代入，得

x  xk
H 3 ( x)  1  2
xk 1  xk

2
 x  xk 1 

 yk
 xk  xk 1 

x  xk 1  x  xk

 1  2
xk  xk 1  xk 1  xk

2

 yk 1

2
 x  xk 1 
 mk
 ( x  xk )
 xk  xk 1 
 x  xk
 ( x  xk 1 )
 xk 1  xk
1
0
0
2

 mk 1

(4.12)
余项 R3 ( x)  f ( x)  H 3 ( x) ，类似（4.5）可得
1 ( 4)
R3 ( x)  f ( )( x  xk ) 2 ( x  xk 1 ) 2 ,   ( xk , xk 1 ).
4!
(4.13)
1
0
1
2.5
2.5.1
分段低次插值
高次插值的病态性质
根据区间 [a, b] 上给出的节点做出的插值多项式 Ln (x),
在次数 n 增加时逼近 f (x) 的精度不一定也增加.
这是因为对任意的插值节点，当 n  时， Ln (x) 不
一定收敛到 f (x) .
考虑函数 f ( x)  1 /(1  x 2 )，它在 [5,5] 上的各阶导数均
存在. 以 [5,5]上的 n  1个等距节点
xk  5  10
k
,
n
(k  0,1, , n)
所构造的拉格朗日插值多项式为
n 1 ( x)
1
.
2
 1 ( x j )
j  0 1  x j ( x  x j )n
n
Ln ( x)  
令 xn 1/ 2  1 ( xn 1  xn ), 则 xn 1/ 2  5  5 ,
2
1
0
3
n
表2-5列出了 n  2,4,,20 时的 Ln ( xn1/ 2 )的计算结果及
在 xn 1/ 2 上的误差 R( xn 1/ 2 ).
表2  5
1
0
4
n
f ( xn 1/ 2 )
Ln ( xn 1/ 2 )
R ( xn 1/ 2 )
2
0.137931
0.759615
 0.621684
4
0.066390
 0.356826
0.423216
6
0.054463
0.607879
 0.553416
8
0.049651
 0.831017
0.880668
10
0.047059
1.578721
 1.531662
12
14
0.045440
0.044334
 2.755000
5.332743
2.800440
 5.288409
16
0.043530
 10.173867
10.217397
18
0.042920
20.123671
 20.080751
20
0.042440
 39.952449
39.994889
2.5
2
1.5
1
0.5
0
1
0
5
-0.5
-5
-4
-3
-2
-1
0
1
2
3
4
5
可见，随 n 的增加，R( xn 1/ 2 ) 的绝对值几乎成倍增加.
这说明当 n   时 Ln在 [5,5]上是不收敛的.
Runge证明了，存在一个常数 c  3.63，使得当 x  c
Ln ( x)  f ( x), 而当 x  c 时 Ln (x) 发散.
时，lim
n 
1
0
6
从图上看到，在 x  5 附近, L10 ( x)与 f ( x)  1 /(1  x 2 )
偏离很远， 这说明用高次插值多项式 Ln (x)近似 f (x) 效
果并不好.
通常不用高次插值，而用分段低次插值.
1
0
7
下图是用Matlab完成的Lagrange插值（附程序）：
1
0
8
附：Lagrange插值程序
1
0
9
n=11; m=61;
x= -5:10/(m-1):5;
y=1./(1+x.^2);
z=0*x;
x0=-5:10/(n-1):5;
y0=1./(1+x0.^2);
y1=lagr1(x0, y0, x);
plot(x, z, ’r’, x, y, ’k:’ ,x, y1, ’r’)
gtext(‘Lagr.’),
gtext(‘y=1/(1+x^2)’)
title(‘Lagrange’)
附：Lagrange插值子程序 lagr1:
1
1
0
function y=lagr1(x0,y0,x)
n=length(x0); m=length(x);
for i=1:m
z=x(i);
s=0.0;
for k=1:n
p=1.0;
for j=1:n
if j~=k
p=p*(z-x0(j))/(x0(k)-x0(j));
end
end
s=p*y0(k)+s;
end
y(i)=s;
end
2.5.2
分段线性插值
由于升高插值多项式的阶数有时并不能达到提高精度
的效果, 所以实际中往往采用分段插值的思想.
分段插值的基本思想是将插值区间划分为若干个小区
间, 然后在每个小区间上做满足一定条件的低阶插值.
所谓分段线性插值就是通过插值点用折线段连接起来
逼近 f (x).
1
1
1
设已知节点 a  x0  x1    xn  b 上的函数值
f 0 , f1 ,, f n , 记 hk  xk 1  xk , h  max hk ,
k
求一折线函数 I h (x) ,
满足：
1. I h ( x)  C[a, b],
2. I h ( xk )  f k
(k  0,1,, n),
3. I h ( x) 在每个小区间 [ xk , xk 1 ]上是线性函数.
则称 I h (x)为分段线性插值函数.
1
1
2
由定义可知 I h (x) 在每个小区间 [ xk , xk 1 ] 上可表示为
x  xk 1
x  xk
I h ( x) 
fk 
f k 1
xk  xk 1
xk 1  xk
( xk  x  xk 1 ),
(5.1)
k  0,1, , n  1.
分段线性插值的误差可利用插值余项（2.15）得到
max
xk  x  xk 1
M2
f ( x)  I h ( x) 
max ( x  xk )( x  xk 1 )
x
2 k  x xk 1
或写成
max f ( x)  I h ( x) 
a  x b
1
1
3
其中 M 2  max f ( x) .
a  x b
M2 2
h ,
8
(5.2)
由此还可以得到
lim I h ( x)  f ( x)
h 0
在 [a, b]上一致成立，故 I h (x) 在 [a, b] 上一致收敛到 f (x).
1
1
4
下图是用Matlab完成的分段线性插值（附程序）：
1
1
5
附：分段线性插值程序
n=11; m=61;
x=-5:10/(m-1):5;
y=1./(1+x.^2);
z=0*x;
x0=-5:10/(n-1):5;
y0=1./(1+x0.^2);
y1=interp1(x0, y0, x);
plot(x, z, ’r’, x, y, ’k:’, x, y1, ’r’)
gtext(‘Piece. –linear.’), gtext(‘y=1/(1+x^2)’)
title(‘Piecewise Linear’)
注：interp1(x0,y0,x)为Matlab中现成的分段线性插值程序.
1
1
6
2.5.3
分段三次埃尔米特插值
分段线性插值函数 I h (x) 的导数是间断的，若在节点
xk (k  0,1,, n)上除已知函数值 f k外还给出导数值
f k  mk (k  0,1,, n), 就可以构造出一个导数连续的分段
插值函数 I h (x) ，满足条件
1. I h ( x)  C1[a, b],
2. I h ( xk )  f k , I h ( xk )  f k (k  0,1,, n),
3. I h ( x)在每个小区间[ xk , xk 1 ] 上是三次多项式.
1
1
7
根据两点三次埃尔米特插值插值多项式（4.12），
I h (x) 在区间 [ xk , xk 1 ] 上的表达式为
 x  xk 1 

I h ( x)  
 xk  xk 1 
2

x  xk
1  2
xk 1  xk


 x  xk
 f k  

 xk 1  xk



2
2

 x  xk 1 
x  xk 1 
 f k 1  
  x  xk  f k
 1  2
xk  xk 1 

 xk  xk 1 
 x  xk
 
 xk 1  xk
2

  x  xk 1  f k1.

上式对于 k  0,1,, n  1 成立.
1
1
8
(5.3)
利用三次埃尔米特插值的余项（4.13），可得误差估计
f ( x)  I h ( x) 
1 4
hk max f ( 4) ( x) ,
384 xk  x xk 1
x  [ xk , xk 1 ],
hk  xk 1  xk , 于是有
定理3
设 f  C 4 [a, b], I h ( x) 为 f (x)在节点
a  x0  x1    xn  b
上的分段三次埃尔米特插值多项式，则有
h4
max f ( x)  I h ( x) 
max f ( 4) ( x) ,
a  x b
384 a xb
其中 h  max ( xk 1  xk ).
1
1
9
0 k  n 1
2.6
三次样条插值
2.6.1
定义3
三次样条函数
若函数 S ( x)  C 2 [a, b], 且在每个小区间 [ x j , x j 1 ]
上是三次多项式，其中 a  x0  x1    xn  b是给定节点，
则称 S (x)是节点 x0 , x1 ,, xn上的三次样条函数.
若在节点 x j 上给定函数值 y j  f ( x j )( j  0,1,, n),
并成立
S(x j )  y j
( j  0,1,, n),
则称 S (x) 为三次样条插值函数.
1
2
1
（6.1）
由于 S (x)在每个小区间 [ x j , x j 1 ] 上有4个待定系数，
共有 n个小区间，所以共有4n 个待定参数.
因为 S ( x )在 [a, b] 上二阶导数连续，所以在节点
x j ( j  1,2, , n  1) 处应满足连续性条件
S ( x j  0)  S ( x j  0),
S ( x j  0)  S ( x j  0).
S ( x j  0)  S ( x j  0),
（6.2）
这些共有 3n  3 个条件，再加上S (x ) 本身还要满足的 n  1
个插值条件，共有 4n  2个条件，还需要2个才能确定 S (x ) .
1
2
2
通常可在区间 [a, b] 端点 a  x0 , b  xn 上各加一个条件
（称为边界条件），常见的边界条件有以下3种：
1. 已知两端的一阶导数值，即
S ( x0 )  f 0,
2.
S ( xn )  f n.
（6.3）
已知两端的二阶导数，即
S ( x0 )  f 0,
S ( xn )  f n,
（6.4）
其特殊情况为
S ( x0 )  S ( xn )  0.
（6.4）‘称为自然边界条件.
1
2
3
（6.4）‘
3. 当 f (x)是以 xn  x0 为周期的周期函数时，则要求
S (x ) 也是周期函数.
这时边界条件应满足
S ( x0  0)  S ( xn  0),
S ( x0  0)  S ( xn  0),
S ( x0  0)  S ( xn  0).
此时插值条件（6.1）中 y0  yn.
这样确定的样条函数 S (x)称为周期样条函数.
1
2
4
（6.5）
2.6.2
样条插值函数的建立
下面利用 S (x ) 的二阶导数值 S ( x j )  M j ( j  0,1,, n)
表示 S (x ) .
由于 S (x ) 在区间 [ x j , x j 1 ] 上是三次多项式，故 S (x)
在 [ x j , x j 1 ] 上是线性函数，可表示为
S ( x)  M j
x j 1  x
hj
 M j 1
x  xj
（6.7）
hj
对 S (x) 积分两次并利用 S ( x j )  y j 及 S ( x j 1 )  y j 1 ，
可定出积分常数，于是得三次样条表达式
1
2
5
S ( x)  M j
( x j 1  x)3
6h j

M j h 2j
 yj 

6

 M j 1
( x  x j )3
6h j
 x j 1  x 
M j 1h 2j

  y j 1 
 h

6
j


( j  0,1,, n  1).
这里 M j ( j  0,1, , n)，是未知的.
1
2
6
 x  xj

 h
j

（6.8）
为了确定 M j ,对 S (x ) 求导得
S ( x)   M j

( x j 1  x)
2
2h j
y j 1  y j
hj

 M j 1
M j 1  M j
6
( x  x j )2
2h j
（6.9）
hj ;
由此可求得
S ( x j  0)  
1
2
7
hj
3
Mj 
hj
6
M j 1 
y j 1  y j
hj
.
类似地可求出 S (x) 在区间 [ x j 1 , x j ] 上的表达式，从而得
S ( x j  0) 
h j 1
6
M j 1 
h j 1
3
Mj 
y j  y j 1
h j 1
,
利用 S ( x j  0)  S ( x j  0) 可得
 j M j 1  2M j   j M j 1  d j
( j  1,2, , n  1),
（6.10）
其中
j 
1
2
8
h j 1
h j 1  h j
,
j 
hj
h j 1  h j
,
j  0,1, , n,
dj  6
f [ x j , x j 1 ]  f [ x j 1 , x j ]
h j 1  h j
 6 f [ x j 1 , x j , x j 1 ].
（6.11）
对第一种边界条件（6.3），可导出两个方程
6

( f [ x0 , x1 ]  f 0), 
h0


6
M n 1  2 M n 
( f n  f [ xn 1 , xn ]).
S ( x0 )  f 0h,n 1 S ( xn )  f n. 
2M 0  M 1 
1
2
9
（6.12）
（6.3）
如果令
0  1,
 n  1,
6
( f [ x0 , x1 ]  f 0),
h0
6
dn 
( f n  f [ xn 1 , xn ]),
hn 1
d0 
那么（6.10）及（6.12）可写成矩阵形式
1
3
0
 2 0
 M 0   d 0 

 M   d 
2

 1 
1
1 
 1

（6.13）

      .
 

 ,2, , n  1),
 j M j 1 2M j 2  j M j1 M
 d j ( dj  1
n 1
n 1 
n 1
n 1 

6 2  M   d 
（6.10）


n
n 
2M 0  M 1 
( f [ x0 , x1 ] 
f 0), n 
h0

（6.12）

6
M n 1  2 M n 
( f n  f [ xn 1 , xn ]).

hn 1

对第二种边界条件（6.4），直接得端点方程
M 0  f 0,
M n  f n.
（6.14）
如果令 0  n  0, d 0  2 f 0, d n  2 f n , 则（6.10）和（6.14）
也可以写成（6.13）的矩阵形式.
S ( x0 )  f 0,
1
3
1
S ( xn )  f n,
（6.4）
对于第三种边界条件（6.5），可得
M 0  M n , n M1  n M n 1  2M n  d n ,
（6.15）
其中
S ( x0  0) hS ( xn  0), S ( x0  0)  S ( xn  0),
n 
0
,
 n  1  n 
S ( x0  0) hSn 1(xhn0  0).
dn  6
hn 1
,
hn 1  h0
f [ x0 , x1 ]  f [ xn 1 , xn ]
,
hn 1  h0
（6.10）和（6.15）可以写成矩阵形式
1
3
2
（6.5）
2

 2



 n
1
2

1  M 1 
2

 n 1

2
n
 d1 
 M   d 
 2 
2 

      .
 


n 1  M n 1   d n 1 
2  M n   d n 
（6.16）
( x j 1  x)3
( x  x j )3
S
( x)  M j
 M j 1
（6.13）和（6.16
）
是关于
M j ( j6
6h
h 0,1, , n) 的三对角
j
j
 j 在力学上解释为细梁在
方程组，M
M j h 2j  x j 1  x  x j 截面处的弯矩，称
M j 1h 2j  x  x j


 yj 
  y j 1 

 h
6  hj M 0   d 0  6

2
j
0
为 S (x
的矩，


方程组（6.13）和（6.16）称为三弯矩方程.
)





 M
2

d 
1
1

( j  0,1,, n  1). 1   1（6.8）
（6.13）
（6.13）和（6.16）的系数矩阵为严格对角占优阵，有


 


 .


唯一解,

 


 n 1 2 n 1  M n 1   d n 1 
求解方法可见第5章第4节追赶法，将解得结果代入
n
2  M n   d n 
（6.8）的表达式即可.
1
3
3
例7
设 f (x)为定义在 [27.7,30] 上的函数，在节点
xi (i  0,1,2,3) 上的值如下：
f ( x0 )  f (27.7)  4.1,
f ( x2 )  f (29)  4.1,
f ( x1 )  f (28)  4.3,
f ( x3 )  f (30)  3.0.
试求三次样条函数 S (x ) ，使它满足边界条件 S (27.7)  3.0,
S (30)  4.0.
1
3
4
解
由（6.11）及（6.12）
h0  0.30,
0  1, 1 
10
,
13
dj  6
h1  h2  1,
2 
1
,
2
1 
d0 
2 
1
,
2
3  1,
6
( f [ x0 , x1 ]  f 0)  46.666,
h0
f [ x j , x j 1 ]  f [ x j 1 , x j ]
h j 1  h j
3
,
13
 6 f [ x j 1 , x j , x j 1 ].
d1  6 f [ x0 , x1 , x2 ]  4.00002, d 2  6 f [ x1 , x2 , x3 ]  （6.11）
2.70000,
6

( f [ x0 , x1 ]  f 0), 
6
h

d3 
( f 3  6f0 [ x2 , x3 ])  17.4.
M n 1 h22 M n 
( f n  f [ xn 1 , xn ]).

hn 1

2M 0  M 1 
由此得矩阵形式的方程组（6.13）为
1
3
5
（6.12）
2
3

13



1
2
10
13
1
2
2
1

M0
   46.666 

  4.00002

M
1

.

1 
   2.7000 
M2
2   M 3  17.4000   17.4000 

2  
求解得
M 0  23.531,
1
3
6
M1  0.395,
M 2  0.830,
M 3  9.115.
代入（6.8）得
13.07278( x  28) 3  14.84322( x  28)  0.21944

3

(
x

27
.
7
)
 14.31358( x  27.7), x  [27.7,28],

3

0
.
06583
(
29

x
)
 4.23417(29  x3 )  0.13833

3
S ( x)  
( x j 1 3 x)
(x  x j )
x  [28,29],
 ( xj  28)  3.96167
S ( x) M
 M j(x1  28),
6h j 3
6h
0.13833(30
 x)  3.96167(30 j x)  1.51917
2
2
 


 x
M
h
x

x
M
h
x],j
3
j 
j 4.51917
j 1
j 1x 
j [ 29
,
30


(
x

29
)
(
x

29
)




 y j 1 
  y j 


 h
6
h
6
j
j




（曲线见图2-6）
( j  0,1,, n  1).
（6.8）
1
3
7
图2-6
1
3
8
例8 给定函数 f ( x) 
1
,
2
1 x
 5  x  5, 节点
xk  5  k (k  0,1,,10), 用三次样条插值求 S10 ( x).
取 S10 ( xk )  f ( xk ) (k  0,1,,10), S10
 (5)  f (5),
 (5)  f (5).
S10
直接上机计算可求出 S10 ( x) 在表2-6所列各点的值.
1
3
9
表2  6
x
1
1 x2
S10 ( x)
L10 ( x)
x
1
1 x2
S10 ( x)
L10 ( x)
 5.0 0.03846 0.03846
0.03846
 2.3 0.15898 0.16115 0.24145
 4.8 0.04160 0.03758
1.80438
 2.0 0.20000 0.20000 0.20000
 4.5 0.04706 0.04248
1.57872
 1.8 0.23585 0.23154 0.18878
 4.3 0.05131 0.04842
0.88808
 1.5 0.30769 0.29744 0.23535
 4.0 0.05882 0.05882 0.05882  1.3 0.37175 0.36133 0.31650
 3.8 0.06477 0.06556  0.20130  1.0 0.50000 0.50000 0.50000
 3.5 0.07547 0.07606  0.22620  0.8 0.60976 0.62420 0.64316
 3.3 0.08410 0.08426  0.10832  0.5 0.80000 0.82051 0.84340
 3.0 0.10000 0.10000 0.10000  0.3 0.91743 0.92754 0.94090
 2.8 0.11312 0.11366
 2.5 0.13793 0.13971
1
4
0
0.19837
0.25376
0
1.0000
1.0000
1.0000
下图是用Matlab完成的样条插值（附程序）：
1
4
1
附：样条插值程序
n=11; m=61;
x=-5:10/(m-1):5;
y=1./(1+x.^2);
z=0*x;
x0=-5:10/(n-1):5;
y0=1./(1+x0.^2);
y1=interp1(x0, y0, x, ’spline’);
plot(x, z, ’r’, x, y, ’k:’, x, y1, ’r’)
gtext(‘Spline’), gtext(‘y=1/(1+x^2)’)
title(‘Spline’)
1
4
2
注：interp1(x0, y0, x, ’spline’)为Matlab中现成的样条插值
程序.
也可以将三种插值结果画在一起：
1
4
3
误差界与收敛性
2.6.3
设 f ( x)  C 4 [a, b], S ( x )为满足第一种或第二
定理4
种边界条件（6.3）或（6.4）的三次样条函数，令
h  max hi , hi  xi 1  xi (i  0,1, , n  1),
0i  n 1
则有估计式
S ( x0 )  f 0,
max f
(k )
( x)  S
S ( x0 )  f 0,
a  x b
S ( xn )  f n.
(k )
（6.3）
( x)  Ck max f
S ( xn )  f n, a  x b
( 4)
( x) h 4  k ,
k  0,1,2,
其中 C0 
1
4
4
5
1
3
, C1 
, C2  .
384
24
8
（6.4）
（6.17）
这个定理不但给出了三次样条插值函数S ( x )的误差
估计.
还说明了当 h  0 时， S ( x )及其一阶导数 S (x)和
二阶导数 S ( x ) 均分别一致收敛于 f ( x )， f (x) 及 f ( x).
1
4
5
1
4
6
1
4
7
1
4
8
1
4
9
1
5
0
1
5
1
高元插值及Matlab实现
1
5
3
1
5
4
1
5
5
 mesh(X,Y,Z)
1
5
6
绘制网格图
>> surf(Z) 绘制表面图
1
5
7
>> surf(X,Y,Z)
1
5
8
1
5
9
1
6
0
1
6
1
单调数据点上的二元插值
1
6
3
1
6
4
1
6
5
1
6
6
1
6
7
1
6
8
1
6
9
1
7
0
1
7
1
1
7
2
1
7
3
1
7
4
1
7
5
1
7
6
1
7
7