数值分析课程报告

插值法和多项式拟合的研究
摘要
在科研和生产实践中，常常需要通过一组测量数据来寻找变量x与y的函数关系近似表达式。

解决这类问题的方法有两种：一种是插值法,另一种是拟合法。

插值法的原理是用一个简单函数逼近被计算函数，然后用该简单函数的函数值近似替代被计算函数的函数值。

拟合法能够是从给定的一组实验数据出发，寻找函数的一个近似表达式,该近似表达式能反映数据的基本趋势而又不一定过全部的点，即曲线拟合。

本文主要介绍拉格朗日插值法、埃尔米特插值法、三次样条插值法以及基于最小二乘法的多项式拟合。

关键词：拉格朗日插值，埃尔米特插值，样条插值，多项式拟合
1方法的意义
在许多实际问题及科学研究中，因素之间往往存在着函数关系，然而，这种关系经常很难有明显的解析表达，通常只是由观察与测试得到一些离散数值。

有时，即使给出了解析表达式，却由于表达式过于复杂，不仅使用不便，而且不易于进行计算与理论分析。

解决这类问题的方法有两种：一种是插值法,另一种是拟合法。

插值法的原理是用一个简单函数逼近被计算函数，然后用该简单函数的函数值近似替代被计算函数的函数值。

它要求给出函数的一个函数表，然后选定一种简单的函数形式，比如多项式、分段线性函数及三角多项式等，通过已知的函数表来确定一个简单的函数()x ϕ作为()f x 的近似，概括地说，就是用简单函数为离散数组建立连续模型。

插值法在实际应用中非常广泛，但是它也有明显的缺陷，一是测量数据常常带有测试误差，而插值多项式又通过所有给出的点，这样就是插值多项式保留了这些误差；二是如果实际得到的数据过多，则必然得到次数较高的插值多项式，这样近似的效果并不理想。

拟合法能够很好的解决这些问题，它从给定的一组实验数据出发，寻找函数的一个近似表达式y=()x ϕ,该近似表达式能反映数据的基本趋势而又不一定过全部的点，即曲线拟合的问题，函数的近似表达式y=()x ϕ称为拟合曲线。

常用最小而二乘法来确定拟合曲线。

2插值法的介绍
2.1 插值法定义
设 f (x )为[a ，b ]上的函数，在互异点n x x x ,...,,10处的函数值分别为 )(),...,(),(10n x f x f x f ,构造一个简单函数 ϕ(x ) 作为函数 f (x ) 的近似表达式y = f (x ) ≈ ϕ(x ),使
)()(i i x f x =ϕ , i ＝0, 1, 2, …,n (1.0) 则称ϕ(x ) 为关于节点n x x x ,...,,10的插值函数；称n x x x ,...,,10 为插值节点；称))((i i x f x , i =1,2,… , n 为插值点；f (x ) 称为被插值函数。

式（1.0）称为插值条件。

这类问题称为插值问题。

插值的任务就是由已知的观测点,为物理量(未知量)建立一个简单的、连续的解析模型，以便能根据该模型推测该物理量在非观测点
处的特性。

常用的插值函数类{()x ϕ}是代数多项式，相应插值问题是代数插值，本文主要介绍三种代数差值：拉格朗日插值，埃尔米特插值和样条插值。

2.2拉格朗日插值
2.2.1两点插值问题
已知)(i i x f y = )1,0(=i ，求满足插值条件i i y x P =)(1 )1,0(=i 的插值多项式)(1x P 。

根据解析几何知识可知，所求的)(1x P 为过点),(),,(1100y x y x 的直线，即：
)()(00
10101x x x x y y y x P ---+= 上式经整理可改写为：
式中1010)(x x x x x l --=，0
101)(x x x x x l --=。

显然{)1(0)0(1)(0===i i x l i ,{)1(1)0(0)(1===i i x l i ,且)(0x l ，)(1x l 为由插值节点唯一确定的线性函数。

)(0x l ，)(1x l 为节点10,x x ，上的一次插值基函数。

可以看出，节点10,x x 上的插值基函数的次数为插值节点个数减一，基函数组中所含的函数个数与插值节点数相同。

而满足i i y x P =)(1 )1,0(=i 的插值多项式)(1x P
就是节点10,x x 上插值基函数的线性组合，其组合系数分别为10,y y 。

这种表示为插值基函数线性组合的一次插值多项式也就是一次拉格朗日插值多项式。

当给定n+1个插值节点后，可类似定义n 次插值基函数，并以此构造n 次拉格朗日多项式。

2.2.2 n 次拉格朗日插值多项式
n 次插值基函数：
)
)...()()...(())...()()...(()(110110n k k k k k k n k k k x x x x x x x x x x x x x x x x x l --------=+-+- )...,1,0(n k = 显然)(x l k 具有以下性质：
性质1，{)
(0)(1)(k i k i x l i k === )...,1,0(n k = 性质2，)(x l k )...,1,0(n k =为由插值节点n x x x ,...,,10唯一确定的n 次函数。

性质3， 基函数组所含的基函数个数与插值节点个数相同。

以n 次插值基函数为基础，可得拉格朗日插值多项式为：
k n
k n k k k k k k n k k n k k k n y x x x x x x x x x x x x x x x x y x l x L ∑∑=+-+-=--------==01101100))...()()...(())...()()...(()()( 00()n n i k k i k i
i k x x y x x ==≠-=-∑∏ 若记∏=+-=n i i n x x x 01)()(ω，则有∏≠=+-='n
k
i i i k k n x x x 01)()(ω。

于是)(x L n 也可写成：
k n k k n k n n y x x x x x L ∑=++'-=011)()()()(ωω
2.2.3拉格朗日插值多项式的余项
)()()(x L x f x R n n -=称为拉格朗日插值多项式的余项，也叫截断误差。

用简单的插值多项式)(x L n 代替复杂的函数)(x f ，这种做法是否有效，取决于截断误差是否满足精度要求。

拉格朗日型余项：
)()!
1()()()()(11x n f x L x f x R n n n n +++=-=ωξ 其中∏=+-=n
i i n x x x 01)()(ω，),(b a ∈ξ且依赖于0x 。

一般情况下，),(b a ∈ξ的具体数值无法知道，但是若能够求出1)1()(max ++≤≤=n n b
x a M x f ，则可以得出插值多项式的截断误差限为：
)(max )!1()(11x n M x R n b
x a n n +≤≤++≤ω 由此看出，)(x R n 的大小除了与1+n M 有关外，还与插值节点有密切关系。

当给定M 个点出的函数值，但仅选用其中)1(1m n n <++个作为插值条件而求某点x 处函数值时，n+1个节点n x x x ,...,,10的选取应该尽可能的接近x ，使得计算的函数值的误差限尽可能的小。

2.3埃尔米特插值
许多实际问题不但要求插值多项式与被插值函数在节点处的函数值相当，而且还需要求其导数值相等。

满足这种要求的插值多项式就是埃尔米特插值多项式。

一般情形的埃尔米特插值问题
一般情形的埃尔米特插值问题是指所满足的插值条件中函数值的个数与导数值的个数相等。

即当函数)(x f 在区间],[b a 上n+1个节点i x ),...,1,0(n i =处的函数值)(i i x f y =及导数值)(i i x f m '=给定时，要求一个次数不超过2n+1的多项式)(12x H n +，使之满足：
),...,1,0()()(1212n i m x H y x H i i n i i n =⎭⎬⎫='=++
这里给出个2n+2个插值条件，可唯一确定一个形式为
12121012...)(++++++=n n n x a x a a x H 的多项式，但是要确定2n+2个系数非常复杂，因此此时可以借用构造拉格朗日插值多项式的基函数。

设)(),(x x j j βα ),...,1,0(n j =为次数不超过2n+1的多项式，且满足：
),...,1,0,()(,0)(0)(,)(n j i x x x x ij i j i j i j ij i j =⎪⎭
⎪
⎬⎫='
=='=δββαδα 则满足上述条件的埃尔米特插值多项式可以写成用插值基函数表示的形式： ∑=++=n j j j i j n m x y x H 012])()([βα
上式满足插值条件。

可以得到基函数)(),(x x j j βα的解析式为：
)(1)(21)(20x l x x x x x j n j k k k j j j ⎥⎥⎥⎦
⎤⎢⎢⎢⎣⎡---=∑≠=α ),...,1,0(n j = )()()(2
x l x x x j j j -=β ),...,1,0(n
j = 式中)(x l j 为拉格朗日插值基函数。

因此埃尔米特插值多项式即为：
j j n j j i n
j j n j k k k j j n m x l x x y x l x x x x x H )()()(1)(21)(2
002012∑∑∑==≠=+-+⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---= 一般情形下的埃尔米特插值多项式的余项为：
)()!
22()()()()(12)22(1212x n f x H x f x R n n n n +++++=-=ωξ 式中),(b a ∈ξ，且与x 有关。

埃尔米特插值的几何意义：曲线)(12x H y n +=与曲线)(x f y =在插值节点处有相同的公共切线。

在带导数的插值问题中，有时插值条件中的函数值个数与导数值个数不等。

这时可以以一般情况的埃尔米特插值多项式为基础，运用待定系数法求出满足插值条件的多项式。

2.4样条插值
在上述方法中，我们根据区间],[b a 上给出的节点得以得到函数)(x f 的插值多项式，但是并非插值多项式的次数越高，逼近函数)(x f 的精度越好，主要原因是因为对于任意的插值节点，当∞→n 时，插值多项式)(x P n 不一定收敛到)(x f 。

这种高次插值不准确的现象称为龙格现象。

为了避免高次插值的缺点，人们常常采用分段插值的方法，即将插值区间分为若干个小区间，在每个小区间上运用前面介绍的插值方法构造低次插值多项式。

采用分段现象插值与分段二次插值，可以构造一个整个连续的函数，而采用分段三次埃尔米特插值则可以构造一个整体上具有一节连续导数的插值函数。

实际问题中，很少给出插值点上的导数值。

三次样条插值就是在只给出插值点上函数值的情况下，构造一个整体上具有二阶连续导数的插值函数。

2.4.1三次样条插值函数的定义
设区间],[b a 上取n+1个节点b x x x x x a n n =<<<<<=+-11210...，若函数)(x S 满足条件：
（1）在整个区间],[b a 上具有二阶连续导数
（2）在每个小区间],[1i i x x - ),...,1,0(n i =上是x 的三次多项式
（3）i i y x S =)( ),...,1,0(n
i = 则称)(x S 为)(x f 的三次样条插值。

)(x S 的边界条件为：
（1）给定两端点处的一阶导数值，记为：
00)(m x S ='，n n m x S =')(
（2）给定两端点处的二阶导数值，记为：
00)(M x S =''，n n M x S ='')(
此外，对0)()(0=''=''n x S x S 的边界条件称为自然边界条件。

2.4.2三次样条插值函数的构造
构造三次样条插值函数，就是要写出它在子区间],[1i i x x - ),...,1,0(n i =上的表达式，记为)(x S i ),...,1,0(n i =。

一、用节点处一阶导数表示的三次样条插值函数
记节点处的一阶导数值为i i m x S =')( ),...,1,0(n i =，若已知i m 后，则)(x S 在],[1i i x x - ),...,1,0(n i =上就是满足条件
11)(--=i i y x S ，i i y x S =)(，),...,1,0(n i =
11)(--='i i m x S ，i i m x S =')( ),...,1,0(n i =
的三次埃尔米特插值多项式。

步骤如下：
（1）根据)(x S ，)(x S '在内阶段的连续性及插值条件，运用],[1i i x x -上的二点三次埃尔米特插值多项式，写出)(x S 用i m ),...,1,0(n i =表示的形式。

（2）利用)(x S ''在内节点i x )1,...,1,0(-=n i 的连续性及边界条件，导出含i m ),...,1,0(n i =的n+1阶线性方程组。

（3）求解含i m ),...,1,0(n i =的线性方程组，将得到的i m 代入],[1i i x x -上的二点三次埃尔米特插值多项式。

具体公式为：
i i i i i i i i i i i i i i i y x x x x x x x x y x x x x x x x x x S 2
111121112121)(⎪⎪⎭⎫ ⎝⎛--⎪⎪⎭⎫ ⎝
⎛--++⎪⎪⎭⎫ ⎝⎛--⎪⎪⎭⎫ ⎝⎛--+=------- +i i i i i i i i i i m x x x x x x m x x x x x x 211121)()(⎪⎪⎭⎫ ⎝⎛---+⎪⎪⎭⎫ ⎝⎛------- i m 的确定由)(x S 得二阶导数在内节点连续的条件得到，即：
i i i i i i f m M m m =+++-112λ )1,...,1,0(-=n
i 该式即为n+1阶线性方程组。

其中：
⎪⎪⎪⎩
⎪⎪⎪⎨⎧-+-=-=+=-=+-++++-)(311111111i i i i i i i i i i i i i i i i i i h y y h y y M f M h h h x x h λλλ )1,...,1,0(-=n i 求解该线性方程组还需两个边界条件。

对于第一种边界条件，将其代入，得到只含n-1个未知数的线性方程组：
⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛--=⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡--------n n n n n n n n n m M f f f m f m m m m M M M 11
2201112211222212222 λλλλ 对于第二种边界条件，将其代入，得到n+1阶线性方程组：
⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----n n n n n n f f f f f m m m m m M M M 1210121011
221121222
12 λλλ 其中：
⎪⎪⎩
⎪⎪⎨⎧+-=--=-n n n n n n M h h y y f M h h y y f 23231011010
二、用节点处二阶导数表示的三次样条插值函数
记节点处的二阶导数值为i i m x S ='')( ),...,1,0(n i =，由于)(x S ''在],[1i i x x - ),...,1,0(n i =上是x 的线性函数，因此构造以节点处二阶导数表示的三次埃尔米特插值多项式步骤如下：
（1）根据)(x S ''在内节点的连续性及为线性函数的特点，将)(x S ''表示为线性函数。

在根据)(x S 在内节点的连续性及插值条件，写出)(x S 用i M ),...,1,0(n i =表示的形式。

（2）利用)(x S '在内节点i x )1,...,1,0(-=n i 的连续性及边界条件，导出
含 i M ),...,1,0(n i =的n+1阶线性方程组。

（3）求解含i M ),...,1,0(n i =的线性方程组，将得到的i M 代入],[1i i x x -上的 )(x S 的表达式，即得到二点三次埃尔米特插值多项式。

具体公式为：
i
i i i i i i i i i i i i i i i i h x x h M y h x x h M y M h x x M h x x x S 122113113)6()6(6)(6)()(-------+--+-+-= i M 的确定由)(x S 得二阶导数在内节点连续的条件得到，即：
i i i i i i f M M M M =+++-112λ )1,...,1,0(-=n
i 该式即为n+1阶线性方程组。

其中：
⎪⎪⎪⎩
⎪⎪⎪⎨⎧---+=-=+=-=+-+++++-)(6111111111i i i i i i i i i i i i i i i i i i h y y h y y h h f M h h h M x x h λ )1,...,1,0(-=n i 解该方程还需要两个边界条件。

对于第一种边界条件，将其代入，得到只含n-1个未知数的线性方程组：
⎪⎪⎪⎪⎪⎪⎭
⎫ ⎝⎛=⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-----n n n n n n f f f f M M M M M M 11012211111212212 λλ 式中：
⎪⎪⎩
⎪⎪⎨⎧--=--=-)(6)(61010110n n n n n n h y y m h f m h y y h f 对于第二种边界条件，将其代入，得到n+1阶线性方程组：
⎪⎪⎪⎪⎪⎪⎭
⎫ ⎝⎛--=⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡--------n n n n n n n n n M f f f M M f M M M M M M M 112201112211222212222λλλλ 对于三次样条插值函数来说，当插值节点逐渐加密时，可以证明：不但样条
插值函数收敛于函数本身，而且其导数也收敛于函数的导数。

3.曲线拟合
3.1 用最小二乘法求解矛盾方程组
工程实际中许多问题都可归结为矛盾方程组，实际中需要寻求矛盾方程组的未知数n x x x ,...,,21的一组取值，它使得偏差的绝对值之和∑=N
i i 1δ尽可能的小，为了
便于分析计算和应用，常采用使偏差的平方和2
1112
)(∑∑∑===-==N i n j i j ij N i i b x a Q δ达到最小值，这一条件称为最小二乘原则。

按最小二乘原则来选择n x x x ,...,,21的一组取值的方法称为求解矛盾方程组的最小二乘法，符合条件的n x x x ,...,,21的一组取值称为矛盾方程组的最小二乘解。

把Q 看出n 个自变量n x x x ,...,,21的二次函数，记为),...,,(21n x x x f Q =，因此，求矛盾方程组的最小二乘解就是二乘函数的最小值点。

若矛盾方程组i
j n
j ij b x a =∑=1的系数矩阵A 的秩为n ，则二次函数),...,,(21n x x x f Q =一定存在最小值。

即矛盾方程组的最小二乘解存在，且正则方程组b A Ax A T T =有唯一解，此解就是矛盾方程组的最小二乘解。

3.2多项式拟合
设通过测量得到函数)(x f y =的一组数据为：)(i i x f y = ),...,1,0(n i =，求一个次数低于N-1的多项式m m x a x a x a a x y ++++==...)(2210ϕ )1(-<N m （其中
m a a a ,...,,10待定）
，使其最好的拟合这组数据，最好的标准时：使得)(x ϕ在i x 的偏差i i i y x -=)(ϕδ的平方和2
112
)]([∑∑==-==N i i i N i i y x Q ϕδ达到最小。

m a a a ,...,,10的解即正则方程组b A Ax A T T =的解。

⎪⎪⎪⎪⎩⎪⎪⎪⎪
⎨⎧=++++=++++=++++∑∑∑∑∑∑∑∑∑∑∑∑∑∑===+=+===+=======N i i m i N i m
i m N i m i N i m i N
i m i N
i i
m i N i m i m N i i N i i N i i N
i i
N i m
i m N i i N i i y x x a x a x a x a y x x a x a x a x a y x a x a x a N a 1121221111
01
11132121101
1122110.........
由该正则方程组求得的唯一解代入拟合多项式)(x y ϕ=，即为所求。

5.数值试验
表1是1971年到1990年我国总人口的统计数字，试根据1971年到1985年这15年人口统计数字用下面几种方法预测未来20年的人口数字，并用图示的方法比较1986年到1990年间预测人口数字与实际统计数字的差异，在你所使用的几种预测方法中找出一组较为合理的预测方法。

（1）指数形式bx ae y =；
（2）拉格朗日插值、埃尔米特插值、样条插值； （3）三次多项式拟合； （4）四次多项式拟合
表1 人口统计数字
年份 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 人口 8.5229 8.7177 8.9211 9.0859 9.2420 9.3717 9.4974 9.6259 9.7542 9.8705 年份 1981
1982
1983
1984
1985
1986
1987
1988
1989
1990
人口
10.0072 10.1654 10.3008 10.4357 10.5851 10.7507 10.9300 11.1026 11.2704 11.4333
数值试验结果：1.指数形式拟合图形为：
2.三次多项式拟合的效果：
3.四次多项式拟合的效果
结论：基于最小二乘法的多项式拟合能取得良好的函数逼近效果。

阶数越高，拟合效果越好。

参考文献
[1] 王正林，龚纯，等．精通MATLAB科学计算．北京：电子工业出版社，2009:135
—157
[2]任玉杰．数值分析及其MATLAB实现［M］．北京：高等教育出版社，2007：
584—642
[3]李信真，车刚明，等．计算方法[M]．西安：西北工业大学出版社，2010：
150—167
[4] 金一庆，陈越．数值方法[M].北京：机械工业出版社，2000：165—176
附录
1.按指数形式拟合的程序：
x=[8.5229,8.7177,8.9211,9.0859,9.2420,9.3717,9.4974,9.6259,9.7542 ,9.8705,10.0072,10.1654,10.3008,10.4357,10.5851,10.7507,10.9300,1 1.1026,11.2704,11.4333];
y=[1971,1972,1973,1974,1975,1976,1977,1978,1979,1980,1981,1982,19 83,1984,1985,1986,1987,1988,1989,1990];
p=polyfit(x,log(y),1);
a=exp(p(2));
b=p(1);
y_n=a*exp(b*x);
plot(x,y,x,y,'*',x,y_n),grid,axis([8.5,11.8,1970,1991]),legend('实际统计曲线','实际统计值','指数曲线')
2.按三次多项式拟合的程序：
x=[8.5229,8.7177,8.9211,9.0859,9.2420,9.3717,9.4974,9.6259,9.7542 ,9.8705,10.0072,10.1654,10.3008,10.4357,10.5851,10.7507,10.9300,1 1.1026,11.2704,11.4333];
y=[1971,1972,1973,1974,1975,1976,1977,1978,1979,1980,1981,1982,19 83,1984,1985,1986,1987,1988,1989,1990];
p=polyfit(x,y,3);
y1=polyval(p,x);
plot(x,y,x,y,'*',x,y1),grid,axis([8.5,11.8,1970,1991]),legend('实际统计曲线','实际统计值','三次多项式拟合')
3.按四次多项式拟合的程序：
x=[8.5229,8.7177,8.9211,9.0859,9.2420,9.3717,9.4974,9.6259,9.7542 ,9.8705,10.0072,10.1654,10.3008,10.4357,10.5851,10.7507,10.9300,1 1.1026,11.2704,11.4333];
y=[1971,1972,1973,1974,1975,1976,1977,1978,1979,1980,1981,1982,19 83,1984,1985,1986,1987,1988,1989,1990];
p=polyfit(x,y,4);
y1=polyval(p,x);
plot(x,y,x,y,'*',x,y1),grid,axis([8.5,11.8,1970,1991]),legend('实际统计曲线','实际统计值','四次多项式拟合')。