第十章数据拟合与插值

第十章数据拟合与插值§10.1引言在解决实际问题的生产（或工程）实践和科学实验过程中，通常需要通过研究某些变量之间的函数关系来帮助我们认识事物的内在规律和本质属性，而这些变量之间的未知函数关系又常常隐含在从试验、观测得到的一组数据之中。因此，能否根据一组试验观测数据找到变量之间相对准确的函数关系就成为解决实际问题的关键。例如在工程实践和科学实验中，常常需要从一组试验观测数据(xi,yj)，i=0,1,.,n之中找到自变量x与因变量y之间的函数关系，一般可用一个近似函数y=f(x)来表示。函数y=f(x)的产生办法因观测数据和要求不同而异，通常可采用数据拟合与函数插值两种办法来实现。数据拟合主要是考虑到观测数据受随机观测误差的影响，进而寻求整体误差最小、能较好反映观测数据的近似函数y=f(x)，此时并不要求所得到的近似函数y=f(x)满足yi=f(xi)，i=0,1,.,n。函数插值则要求近似函数y=f(x)在每一个观测点xi处一定要满足yi=f(xi)，i=0,1,.,n。在这种情况下，通常要求观测数据相对比较准确，即不考虑观测误差的影响。在实际问题中，通过观测数据能否正确揭示某些变量之间的关系，进而正确认识事物的内在规律与本质属性，往往取决于两方面因素。其一是观测数据的准确性或准确程度，这是因为在获取观测数据的过程中一般存在随机测量误差，导致所讨论的变量成为随机变量。其二是对观测数据处理方法的选择，即到底是采用插值方法还是用拟合方法，插值方法之中、拟合方法之中又选用哪一种插值或拟合技巧来处理观测数据。插值问题忽略了观测误差的影响，而拟合问题则考虑了观测误差的影响。但由于观测数据客观上总是存在观测误差，而拟合函数大多数情况下是通过经验公式获得的，因此要正确揭示事物的内在规律，往往需要对大量的观测数据进行分析，尤为重要的是进行统计分析。统计分析的方法有许多，如方差分析、回归分析等。数据拟合虽然较有效地克服了随机观测误差的影响，但从数理统计的角度看，根据一个样本计算出来的拟合函数（系数），只是拟合问题的一个点估计，还不能完全说明其整体性质。因此，还应该对拟合函数作区间估计或假设检验，如果置信区间太大或包含零点，则由计算得到的拟合函数系数的估计值就毫无意义。这里所采用的统计分析方法就是所谓的回归分析。另外还可用方差分析的方法对模型的误差作定量分析。对于插值方法，本章在第二节中简单介绍最常用的插值法的基本结论及其Matlab实现问题。由于数据拟合问题必须作区间估计或假设检验，所以除了在本章第三节、第四节中介绍最基本的数据拟合方法——最小二乘法的基本结论及其Matlab实现问题外，我们在第五节中专门介绍了对数值拟合问题进行区间估计或假设检验的统计方法，即介绍回归分析方法及其Matlab实现。数据处理问题通常情况下只是某个复杂实际问题的一个方面或部分内容，因而这里所介绍的数据处理方法——函数插值和数据拟合的方法（包括回归分析）通常只能解决实际问题中的部分问题——计算问题。一般来说，对实际问题进行数学建模需要用到多方面知识，只有很少的情况下可以单独使用本章所介绍的内容，故我们只在本章最后一节以修改后的美国91年数学建模A题为例说明如何使用数值计算知识建立数学模型，从而解决实际问题的方法。§10.2插值方法在实际问题中所遇到的插值问题一般分为一维插值问题和二维插值问题。本节主要介绍最常用的一维插值方法及其一些简单结果。一维插值问题的数学描述为：已知某一函数y=g(x)（g(x)的解析表达式可能十分复杂，也可以是未知的）在区间[a,b]上n+1个互异点xj处的函数值yj，j=0,1,.,n，还知道g(x)在[a,b]上有若干阶导数，如何求出g(x)在[a,b]上任一点x的近似值。一维插值方法的基本思想是：根据g(x)在区间[a,b]上n+1个互异点xj（称为节点）的函数值yj，j=0,1,.,n，求一个足够光滑、简单便于计算的函数f(x)（称为插值函数）作为g(x)的近似表达式，使得f(xj)=yj，j=0，1，…，n。（10.1）然后计算f(x)在区间[a,b](称为插值区间)上点x（称为插值点）的值作为原函数g(x)（称为被插函数）在此点的近似值。求插值函数f(x)的方法称为插值方法，（10.1）称为插值条件。代数多项式比较简单，常用多项式作为插值函数。一插值多项式的存在唯一假设f(x)是一个满足插值条件（10.1）的次数不超过n的代数多项式，即f(x)=a0+a1x+.+anxn（10.2）为满足（10.1）的插值函数，则f(x)的n+1个待定系数a0，a1，…，an满足ìa+ax+ax2+.+axn=y,01020n00.2n.a0+a1x1+a2x1+.+anx1=y1,.....2n.a0+a1xn+a2xn+.+anxn=yn.（10.3）记此方程组的系数矩阵为A，则1xx2.xn0001xx2.xn111det(A)=.......1xx2.xnnnn是范德蒙(Vandermonde)行列式。当x0,x1,.,xn互不相同时，此行列式值不为零。因此方程组（10.3）有唯一解。这表明，只要n+1个插值节点x0,x1,.,xn互异，满足插值条件（10.1）的插值多项式（10.2）存在唯一。从几何上看，上述多项式插值就是过n+1个数据点(xj,yj)，j=0，1，…，n作一条多项式曲线y=f(x)近似原有曲线y=g(x)。当x.[a,b]且x1xj(j=0,1,.,n)时，g(x).f(x)，称被插函数g(x)与插值函数（多项式）f(x)之间的差Rn(x)=g(x)-f(x)为插值函数（多项式）f(x)的截断误差，或插值余项。显然，插值多项式的存在唯一与节点01x,x,.,xn的次序无关。二拉格朗日插值法1拉格朗日插值的基本结果用多项式函数（10.2）作为插值函数时，希望通过解方程组（10.3）而得到待定系数a0，na1，…，a的做法当n比较大时是不现实的。方便的方法是先构造一组基函数li(x)=(x(ix--xx0)(0))..((xxi--xxii--11)()(xxi--xxi+i+11))..((xx-i-xnx)n)，i=0，1，…，n。显然li(x)是n次多项式且满足li(xj)=í.ì10,,jj=1ii,,i,j=0,1,.,n。称n次多项式li(x)为节点x0,x1,.,xn上的n次插值基函数。令Ln(x)=.(n)yili(x)i=0，则多项式Ln(x)显然满足插值条件Ln(xi)=yi(i=0,1,.,n)，从而Ln(x)为插值多项式。插值多项式Ln(x)又称为n次拉格朗日（Lagrange）插值多项式，由方程组（10.3）解的存在唯一性，n+1个互异节点的n次拉格朗日插值多项式Ln(x)存在唯一。当g(x)在[a,b]上充分光滑时，利用罗尔（Rolle）定理可推出：对于任意x.[a,b]，插值多项式Ln(x)的余项(n+1)g(x)nnn+1R(x)=g(x)-L(x)=(n+1)!w(x),x.(a,b)，wn+1(x)=.(n)(x-xj)其中j=0。若可以估计(n+1)|g(x)|￡Mn+1则可以得到n次拉格朗日插值的余项估计(，)nn+1|R(x)|￡(nM+n+11)!|w(x)|。实际上，因为Mn+1常难以确定，所以上式并不能给出精确的误差估计。2拉格朗日插值的Matlab实践Matlab中没有现成的拉格朗日插值函数，必须编写一个M文件实现拉格朗日插值。设n个节点数据以数组x0,y0输入（注意Matlab的数组下标从1开始），m个插值点以数组x输入，输出数组y为m个插值。编写一个名为lagrange.m的M文件：functiony=lagrange(x0,y0,x);n=length(x0);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;forj=1:nifj~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;end三分段线性插值法1高次插值多项式的龙格（Runge）现象用拉格朗日插值多项式Ln(x)作为区间[a,b]上连续函数g(x)的近似函数，在大多数情况下，Ln(x)的次数越高，逼近g(x)的效果就越好，即误差|Rn(x)|越小。但是对于高次多项式插值问题而言，往往会造成插值多项式Ln(x)的收敛性与稳定性变差，反而逼近效果不理想，甚至发生龙格现象，即在理论上并不能保证当n趋于无穷大时Ln(x)在[a,b]上处处收敛于g(x)。关于这一点，在20世纪初就为龙格（Runge）所发现：在[-1,1]上用n+1个等距节点作被插值函数y(x)=1(1+25x2)的插值多项式Ln(x)，则随着n的增大，Ln(x)振荡越来越大。计算结果与理论证明表明，当n趋于无穷大时，Ln(x)在区间中部收敛于y(x)，但是对于满足条件0.726…￡|x|1的x，Ln(x)并不收敛于y(x)。上述例子说明，使用高次多项式插值是危险的，因此在实际计算中不能使用高次插值。这反而启发我们使用分段插值方法，即将区间[a,b]分成一些小区间，在每一个小区间上用低次多项式进行插值，在整个插值区间[a,b]上就得到一个分段低次多项式插值函数。区间的剖份可以是任意的，各小区间上插值多项式的次数的选取也可按具体问题的要求而选择。分段低次多项式插值通常有较好的收敛性和稳定性，算法简单，但插值函数光滑性变差。常用的分段多项式插值法有两类：一类是下面将要介绍的分段线性插值法；另一类是后面将要介绍的三次样条插值法。2分段线性插值法假设区间[a,b]上的连续函数g(x)在n+1个节点a=x0x1.xn=b上的函数值g(xj)=yjj=0，1，…，n。则得到xy平面上的n+1个数据点(xj,yj)。连接相邻数据点(xj-1,yj-1)、((，)xj,yj)得到n条线段，它们组成一条折线。把区间[a,b]上这n条折线段表示的函数称为被插函数g(x)关于这n+1个数据点的分段线性插值函数，记作I(x)，则I(x)具有如下性质：(1)I(x)可以分段表示，在每个小区间[xj-1,xj]上，它是线性函数，即x-xjx-xj-1I(x)=yj-1xj-1-xj+yjxj-xj-1，xj-1￡x￡xj；jj(2)I(x)=yj=0，1，…，n；(3)I(x)在[a,b上连续。若构造插值基函数](，)ìx-x.xi-xii--11,x.[xi-1,xi](i=0时舍去）,.li(x)=.íxxi--xxii++11，x.[xi,xi+1]（i=n时舍去）,..0，其它,..I(x)=yl(x)j=0则.(n)jj。当g(x)在[a,b]上连续时，分段线性插值函数I(x)具有良好的收敛性，即limI(x)=g(x)h.0，x.[a,b]，而且当g(x)在[a,b]上二阶导数连续时，对于任意x.[a,b]有|R(x)|=|g(x)-I(x)|￡M82h2，h=max{xj-xj-1}M2=max{|g''(x)|}其中1￡j￡n，a￡x￡b。用I(x)计算x点的插值时，只用到x左右的两个节点，计算量与节点个数n+1无关。但n越大，分段越多，插值误差越小。实际上用数据点作插值计算时，分段线性插值就足够了，如数学、物理中用的特殊函数表，数理统计中用的概率分布表等。3分段线性插值的Matlab实现用Matlab实现分段线性插值不需要编制函数程序，Matlab中有现成的一维插值函数interp1。y=interp1(x0,y0,x,'method')method指定插值的方法，默认为线性插值。其值可为：'nearest'最近项插值'linear'线性插值'spline'立方样条插值'cubic'立方插值。所有的插值方法要求x0是单调的。当x0为等距时可以用快速插值法，使用快速插值法的格式为'*nearest'、'*linear'、'*spline'、'*cubic'。四三次样条插值法1三次样条的基本结果分段线性插值函数在节点处的一阶导数一