数据拟合与插值
第十章第十章 数据拟合与插值数据拟合与插值 §§10.1 引言引言 在解决实际问题的生产（或工程）实践和科学实验过程中，通常需要通过研究某些变量之 间的函数关系来帮助我们认识事物的内在规律和本质属性，而这些变量之间的未知函数关系又 常常隐含在从试验、 观测得到的一组数据之中。 因此，能否根据一组试验观测数据找到变量之间 相对准确的函数关系就成为解决实际问题的关键。 例如在工程实践和科学实验中，常常需要从一组试验观测数据 ),( ji yx ， ni, 1 , 0= 之 中找到自变量x与因变量 y 之间的函数关系，一般可用一个近似函数)(xfy =来表示。函数 )(xfy = 的产生办法因观测数据和要求不同而异，通常可采用数据拟合与函数插值两种办法 来实现。 数据拟合主要是考虑到观测数据受随机观测误差的影响，进而寻求整体误差最小、能较好 反映观测数据的近似函数 )(xfy = ，此时并不要求所得到的近似函数 )(xfy = 满足 )( ii xfy = ， ni, 1 , 0= 。 函数插值则要求近似函数)(xfy =在每一个观测点 i x处一定要满足)( ii xfy = ， ni, 1 , 0= 。在这种情况下，通常要求观测数据相对比较准确，即不考虑观测误差的影响。 在实际问题中，通过观测数据能否正确揭示某些变量之间的关系，进而正确认识事物的内 在规律与本质属性，往往取决于两方面因素。其一是观测数据的准确性或准确程度，这是因为 在获取观测数据的过程中一般存在随机测量误差，导致所讨论的变量成为随机变量。其二是对 观测数据处理方法的选择，即到底是采用插值方法还是用拟合方法，插值方法之中、拟合方法 之中又选用哪一种插值或拟合技巧来处理观测数据。插值问题忽略了观测误差的影响，而拟合 问题则考虑了观测误差的影响。但由于观测数据客观上总是存在观测误差，而拟合函数大多数 情况下是通过经验公式获得的，因此要正确揭示事物的内在规律，往往需要对大量的观测数据 进行分析，尤为重要的是进行统计分析。 统计分析的方法有许多，如方差分析、 回归分析等。 数 据拟合虽然较有效地克服了随机观测误差的影响，但从数理统计的角度看，根据一个样本计算 出来的拟合函数（系数），只是拟合问题的一个点估计，还不能完全说明其整体性质。因此， 还应该对拟合函数作区间估计或假设检验，如果置信区间太大或包含零点，则由计算得到的拟 合函数系数的估计值就毫无意义。 这里所采用的统计分析方法就是所谓的回归分析。 另外还可用 方差分析的方法对模型的误差作定量分析。 对于插值方法，本章在第二节中简单介绍最常用的插值法的基本结论及其Matlab实现问题。 由于数据拟合问题必须作区间估计或假设检验，所以除了在本章第三节、第四节中介绍最基本 的数据拟合方法最小二乘法的基本结论及其Matlab实现问题外，我们在第五节中专门介绍 了对数值拟合问题进行区间估计或假设检验的统计方法，即介绍回归分析方法及其Matlab实现。 数据处理问题通常情况下只是某个复杂实际问题的一个方面或部分内容，因而这里所介绍 的数据处理方法函数插值和数据拟合的方法（包括回归分析）通常只能解决实际问题中的 部分问题计算问题。一般来说，对实际问题进行数学建模需要用到多方面知识，只有很少 的情况下可以单独使用本章所介绍的内容，故我们只在本章最后一节以修改后的美国91年数学 建模A题为例说明如何使用数值计算知识建立数学模型，从而解决实际问题的方法。 §§10.2 插值方法插值方法 在实际问题中所遇到的插值问题一般分为一维插值问题和二维插值问题。本节主要介绍最 常用的一维插值方法及其一些简单结果。 一维插值问题的数学描述为：已知某一函数 )(xgy = （ )(xg 的解析表达式可能十分复杂， 也可以是未知的）在区间,ba上 1+n 个互异点 j x 处的函数值 j y ， nj, 1 , 0= ，还知道 )(xg 在 ,ba 上有若干阶导数，如何求出 )(xg 在 ,ba 上任一点x的近似值。 一维插值方法的基本思想是：根据)(xg在区间, ba上 1+n 个互异点 j x （称为节点）的 函数值 j y ， nj, 1 , 0= ，求一个足够光滑、简单便于计算的函数 )(xf （称为插值函数）作 为 )(xg 的近似表达式，使得 jj yxf=)( ， 0=j ，1，n。 （10.1） 然后计算)(xf在区间 ,ba (称为插值区间)上点x（称为插值点）的值作为原函数 )(xg （称 为被插函数）在此点的近似值。求插值函数 )(xf 的方法称为插值方法，（10.1）称为插值条件。 代数多项式比较简单，常用多项式作为插值函数。 一 插值多项式的存在唯一 假设 )(xf 是一个满足插值条件（10.1）的次数不超过n的代数多项式，即 n nx axaaxf+= 10 )( （10.2） 为满足（10.1）的插值函数，则 )(xf 的 1+n 个待定系数 0 a ， 1 a ， n a 满足 =+ =+ =+ . , , 2 210 11 2 12110 00 2 02010 n n nnnn n n n n yxaxaxaa yxaxaxaa yxaxaxaa （10.3） 记此方程组的系数矩阵为A，则 n nnn n n xxx xxx xxx A 2 1 2 11 0 2 00 1 1 1 )det(= 是范德蒙(Vandermonde)行列式。当 n xxx, 10 互不相同时，此行列式值不为零。因此方程组 （10.3）有唯一解。 这表明，只要 1+n 个插值节点 n xxx, 10 互异，满足插值条件（10.1）的 插值多项式（10.2）存在唯一。 从几何上看，上述多项式插值就是过 1+n 个数据点 ),( jj yx ， j =0，1，n作一条多 项式曲线 )(xfy = 近似原有曲线 )(xgy = 。 当 ,bax 且 j xx ), 1 , 0(nj= 时， )(xg)(xf ，称被插函数 )(xg 与插值函数 （多项式） )(xf 之间的差 )()()(xfxgxRn= 为插值函数（多项式） )(xf 的截断误差，或插值余项。显然，插值多项式的存在唯一与节点 n xxx, 10 的次序无关。 二 拉格朗日插值法 1 拉格朗日插值的基本结果 用多项式函数（10.2）作为插值函数时，希望通过解方程组（10.3）而得到待定系数 0 a ， 1 a ， n a 的做法当n比较大时是不现实的。方便的方法是先构造一组基函数 )()()( )()()( )( 110 110 niiiiii nii i xxxxxxxx xxxxxxxx xl = + + ，i=0，1，n。 显然 )(xli 是n次多项式且满足 = = , 1 , 0 )( ij ij xl ji nji, 1 , 0,= 。 称n次多项式 )(xli 为节点 n xxx, 10 上的n次插值基函数。令 = = n i iin xlyxL 0 )()( ， 则多项式 )(xLn 显然满足插值条件 ), 1 , 0()(niyxL iin = ，从而 )(xLn 为插值多项式。 插 值多项式 )(xLn 又称为n次拉格朗日（Lagrange）插值多项式，由方程组（10.3）解的存在唯 一性， 1+n 个互异节点的n次拉格朗日插值多项式 )(xLn 存在唯一。 当 )(xg 在 ,ba 上充分光滑时，利用罗尔（Rolle）定理可推出：对于任意 ,bax ， 插值多项式 )(xLn 的余项 ),(),( )!1( )( )()()( 1 )1( bax n g xLxgxR n n nn + = + + ， 其中 = + = n j jn xxx 0 1 )()( 。若可以估计 1 )1( | )(| + + n n Mg ， 则可以得到n次拉格朗日插值的余项估计 | )(| )!1( | )(| 1 1 x n M xR n n n+ + + 。 实际上，因为 1+n M 常难以确定，所以上式并不能给出精确的误差估计。 2 拉格朗日插值的Matlab实践 Matlab中没有现成的拉格朗日插值函数，必须编写一个M文件实现拉格朗日插值。 设n个节点数据以数组 0, 0 yx 输入（注意Matlab的数组下标从1开始），m个插值点以数 组x输入，输出数组 y 为m个插值。编写一个名为lagrange.m的M文件： function y=lagrange(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 三 分段线性插值法 1 高次插值多项式的龙格（Runge）现象 用拉格朗日插值多项式 )(xLn 作为区间 ,ba 上连续函数 )(xg 的近似函数，在大多数情况 下， )(xLn 的次数越高，逼近 )(xg 的效果就越好，即误差 | )(|xRn 越小。但是对于高次多项式 插值问题而言，往往会造成插值多项式 )(xLn 的收敛性与稳定性变差，反而逼近效果不理想， 甚至发生龙格现象，即在理论上并不能保证当n趋于无穷大时 )(xLn 在 ,ba 上处处收敛于 )(xg 。关于这一点，在20世纪初就为龙格（Runge）所发现：在 1 , 1 上用 1+n 个等距节点作 被插值函数 )251 (1)( 2 xxy+= 的插值多项式 )(xLn ，则随着n的增大， )(xLn 振荡越来越大。 计算结果与理论证明表明，当n趋于无穷大时， )(xLn 在区间中部收敛于 )(xy ，但是对于满 足条件0.726 | x nFF 时，拒绝 0 H ；否则就接受 0 H 。 （b）t检验法 当 0 H 成立时，由T分布的定义有 e n i i xx T )( 1 1 2 = = )2( nt 。 因此，对于给定的显著水平，用T统计量检验 0 H ，有置信水平为 1 的临界值 )2, 1 ( 2/1 nt ,从而t检验法的检验准则为： 当 )2(| 2/1 ntT 时，拒绝 0 H ；否则就接受 0 H 。 2回归系数的区间估计 令 = = n i i xxV 1 2 )( ，则由 )/1() ( 2 00 Vxn e + )2( nt ， )() ( 11 V e )2( nt 得到在置信水平 1 下， 0 和 1 的置信区间分别为 Vxnnt e /1)2( 2 2/10 + ， Vxnnt e /1)2( 2 2/10 + ， Vnt e/ )2( 2/11 ， Vnt e/ )2( 2/11 + 。 再由 2 / e Q )2( 2 n ，还可得在置信水平 1 下， 2 的置信区间为 )2( 2 2/1 n Qe ， )2( 2 2/ n Qe 。 3） 预测与控制 当检验结果拒绝了 0: 10 =H ，接下来的问题是如何利用回归方程 xY 10 += 进行预 测和控制。预测就是对固定的x值预测相应的 y 值，控制就是通过控
