数学建模上课.doc

第六部分 插值与拟合实际生活中经常会遇到这样的问题：给定一批数据点，需要确定满足特殊要求的曲线或曲面。如果要求所求曲线（面）通过所给所有数据点，这就是插值问题，在数据教少的情况下，这样做能取得较好的效果。但是，如果数据较多，那么插值函数是一个次数很高的函数，比较复杂，同时，给定的数据一般是由观察测量所得，往往带有随机误差，因而，要求曲线（面）通过所有数据点就既不现实也不必要。如果不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势，可得到更简单实用的近似函数，这就是数据拟合。函数插值与数据拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者在数学方法上是完全不同的。一、 插值问题1、拉格朗日插值若知道函数在互异的两个点和处的函数值和，而想估计该函数在另一点处的函数值，最自然的想法是作过点和点的直线，用作为准确值的近似值，如果认为误差太大，还可增加一点的函数值，即已知在互异的三个点和处的函数值和，可以构造一个过这三点的二次曲线，用作为准确值的近似值。【定义1】一般的，若已知在互异的个点处的函数值，则可以考虑构造一个过这个点的次数不超过的多项式，使其满足 （1）然后用作为准确值的近似值。此方法就叫做插值法，这样构造出来的多项式称为的次拉格朗日插值多项式或插值函数。称点为插值结点，称式(1)为插值条件，含的最小区间叫做插值区间。【定理1】满足插值条件(1)的次数不超过的多项式是存在的而且是唯一的。1.1 线性插值公式已知函数在互异的两个点和处的函数值和，欲求一个次数不超过1的多项式，使其满足：， （2）根据定理1，是存在而且唯一的，称为线性插值函数或一次插值多项式。用点斜式可以写出过点和点的直线方程：，因此，将它写成对称式为 （3）我们称(3)式为拉格朗日线性插值函数或一次拉格朗日插值公式。若引入记号：，则(3)式可以写成： (4)其中满足：我们称为线性插值或一次拉格朗日插值的基函数。例1 根据下表给出的平方根的值，用线性插值计算，149161234解：取最接近的两点为插值节点，运用插值公式（3），得1.2抛物线插值公式 已知函数在三个互异点和处的函数值和，欲求一个次数不超过2的多项式，使其满足： (5) 根据定理1，是存在而且唯一的。下面仿照线性插值时构造插值函数的方法，用基函数的线性组合作出满足(5)的二次插值多项式，此时有三个基函数：,它们都是二次函数。分别满足下列条件： (6) 由上述条件可以推出的表达式：例如：由于是的零点，故必为形式，为待定系数，又由于，代入此式解出，从而得到的表达式为：同理可得：根据以上讨论，我们得到二次插值多项式为： (7)称为抛物线插值函数或二次插值多项式。例2 根据例1的数据，用抛物线法计算的近似值。解：选择与最接近的三点为插值节点，根据抛物线插值公式(7)，有练习：已知在处的函数值，求的近似值。1.3一般情形 一般的次插值问题是构造满足条件(1)的次数不超过的多项式。与构造线性和二次插值多项式类似，次拉格朗日插值公式可表成次插值基函数的线性组合， (8)其中是次多项式，且满足 (9)与前面推导类似，可以由(9)式得到的具体表达式： , (10)为便于书写，引进记号：,取在处的导数，得于是拉格朗日插值公式可写为：1.4插值多项式的余项本段讨论用次插值多项式来近似函数时的误差。记并称为次插值多项式的截断误差，或称为插值余项。【定理2】假设函数在上有阶导数，且 (11)是经过数据点的次数不高于的多项式。则对区间上的任何都存在(依赖于)，使得 (12)其中 (13)插值公式的余项可以给出多项式逼近的误差估计。2、分段插值多项式历来都被认为是最好的逼近工具之一。用多项式作插值函数，就是前面的代数插值。一般情况下，似乎可以靠增加插值结点的数目来改善插值的精度，但插值多项式的次数会随着结点个数的增加而升高，可能造成插值函数的收敛性和稳定性变差，逼近的效果往往是不理想的，甚至会发生龙格振荡现象。（龙格在本世纪初发现：在区间上用个等距结点作插值多项式，使得它在结点的值与函数在对应结点的值相等，当时，插值多项式在区间的中部趋于，但是对于满足条件的，并不趋于在对应点的值。这种现象就叫做龙格现象。）若插值的范围较小（在某个局部），用低次插值往往就能奏效，例如对 在每个子段上用线性插值，即用连接相邻结点的折线逼近所考察的曲线，就能保证一定的逼近效果。这种增加结点，用分段低次多项式插值的化整为零的处理办法称为分段插值法，也就是说不是去寻求整个插值区间上的一个高次多项式，而是把插值区间划分为若干个小区间，在每个小区间上用低次多项式进行插值，在整个插值区间上就得到一个分段插值函数。区间的划分是任意的，各个区间上插值多项式的次数的选取也可按具体问题选择。分段插值法通常有较好的收敛性和稳定性，算法简单，克服了龙格现象，但插值函数不如拉格朗日插值多项式光滑。这类插值大致可分为两类：一类是下面要介绍的局部化的简单分段插值；另 一类是非局部化光滑性较好的分段插值，即后面要介绍的样条插值。 2.1分段线性插值在分段插值中，用得较多的是分段线性插值。【定义1】设在区间上取个结点： (14) 在区间上有二阶导数的函数在上列结点的值为(15)于是得到个数据点。联结相邻两点得条线段，它们组成一条折线，把区间上这条折线表示的函数称为函数关于这个数据点的分段插值函数，记为它有如下性质：(1) 可以用分段函数表示，在区间上连续。(2) 在第段区间上的表达式为 (16)由此构造插值基函数：则 (17)对分段线性插值的余项估计有下列结果：【定理1】设给定结点为在上存在，则对任意的，有，其中2.2三次样条插值 分段线性插值函数在结点的一阶导数一般不存在，光滑性不高，这就导致了样条插值的提出。 在机械制造，航海等工业中，经常有这样的问题：已知一些数据点如何通过这些数据点作一条比较光滑（如二阶导数连续）的曲线呢？ 解决这一问题的方法是：首先把数据点描绘在平面上，再把一根富有弹性的细直条（称为样条）弯曲，使其一边通过这些数据点，用压铁固定细直条的形状，沿样条边绘出一条光滑的曲线。往往要用几根样条，分段完成上述工作，这时应当让连接点也保持光滑。对这一用样条绘出的曲线，进行数学模拟，这样就导出了样条函数的概念。【定义2】设在区间上，已给个互不相同的结点 (18)而函数在这些结点的值。如果分段表示的函数满足下列条件， (1) 在子区间的表达式都是次数不高于3的多项式; (2) ;(3) 在整个区间上有连续的二阶导数。就称为在基点的三次样条插值函数，简称三次样条。 在区间上，是不超过三次的多项式，故是线性函数，令 (19)由拉格朗日线性插值公式，得 (20)其中, 由(20)得 (21)由于满足插值原则，并且在全区间有直到二阶的连续导数，又由(21)知在上存在且为常数，所以由泰勒公式有： (22)将代入(22)，可解得 (23)将(23)代入(22)，得到在上的表达式： (24)由此可见，只要能求出,就能完全确定。下面利用存在且连续的条件推导应满足的关系式。在上，的表达式为由上式利用导数的定义容易计算出在的左导数为：(25)令，整理后得到： (26)再令 (27) (28)则(26)式可写为 (29)即得到含有个未知量的个方程的方程组。要确定未知量还需补充两个条件,它们通常在区间的两端给出，称为边界条件。常用的边界条件有三种：(1) 边值条件：给定端点处的一阶导数，即(2) 边值条件：给定端点处的二阶导数，即已知。特别取时称为自然边界条件，求得的称为自然样条函数。(3) 周期边值条件：当是以为周期的周期函数时，要求也是周期函数，故端点要满足。给出任一种边界条件都可以得到两个独立的方程，将它们与(29)联立就得到个方程，解出之后代入(24)便得到三次样条函数在各个子区间上的表达式，以第二种边界条件为例，为已知数，于是(29)化为： (30)三对角方程组(30)的系数矩阵按行严格对角占优，所以它的行列式不等于零，从而方程组的解存在并且唯一，可用追赶法求解。（追赶法的程序）三次样条函数的计算步骤如下：(1) 由(27),(28)计算;(2) 由方程组(29)结合给定的边界条件求出确定的方程组，并求解;(3) 将代入(24)，得到的分段表达式。以上三种插值方法都是一维插值，它们有如下特点：拉格朗日插值（高次多项式插值）其插值函数在整个区间上是一个解析表达式，便于再次开发利用；曲线光滑；误差估计有表达式；收敛性不能保证（振荡现象）。用于理论分析，实际意义不大。分段线性和三次样条插值（低次多项式插值）：曲线不光滑（三次样条插值已大有改进）；误差估计较难（对三次样条插值）；收敛性有保证。简单实用，应用广泛。除此之外，还有Hermit插值，不少实际问题中不但要求在节点上函数值相等，而且要求导数值也相等，甚至要求高阶导数值也相等，满足这一要求的插值多项式就是Hermite插值多项式。下面只讨论函数值与一阶导数值个数相等且已知的情况。已知个插值节点及其对应的函数值和一阶导数值，则计算插值区域内任意的函数值的Hermite插值公式:其中 二. 用MATLAB解插值问题1、 一维插值MATLAB在一维插值函数interp1中，提供了四种插值方法选择：线性插值，三次样条插值，三次插值和最近邻点插值（linear,spline,cubic,nearest）。Interp1的基本格式为：interp1(x,y,cx, method) 对一组节点进行插值，计算插值点的函数值 其中分别表示为节点向量值和对应的节点函数值，如果为矩阵，则插值对的每一列进行，若的维数超出或的维数，则返回。而method为可选参数，对应于上述四种方法，可从以下四个值中任选一个： nearest-最近邻点插值linear-线性插值spline-三次样条插值cubic-三次插值linear是缺省值。当缺省时按线性插值处理。所有的插值方法要求是单调的，也可能并非连续等距的。当为等距且连续时，可以用快速插值法。例1 正弦曲线的插值示例解：输入命令： x=0:0.1:10; y=sin(x);xi=0:.25:10;yi=interp1(x,y,xi);plot(x,y,o,xi,yi)例2在12h内，每隔1h测量一次温度，温度依次为：5,8,9,15,25,31,30,22,25,27,24.试估计在3.2,6.5,7.1,11.7h时的温度值。解：输入命令： hours = 1:12; temps = 5 8 9 15 25 29 31 30 22 25 27 24; t = interp1(hours,temps,3.2 6.5 7.1 11.7 % 缺省时按线性插值 T = interp1（hours,temps,3.2 6.5 7.1 11.7,spline） % 三次样条插值计算结果： t = 10.2000 30.0000 30.9000 24.9000 T= 9.6734 30.0427 31.1755 25.3820比较发现，样条插值和线性插值的结果不同。因为插值是一个估计或猜测的过程，应用不同的估计规则将导致不同的结果。一个最常用的样条插值是对数据进行平滑。即给定一组数据，使用样条插值在更细的间隔内求值。例3 在例3的条件下，每隔1/10h估计一次温度值。解：由于数据点太多，最好用图形表示。输入命令： hours = 1:12; temps = 5 8 9 15 25 29 31 30 22 25 27 24; h =1:0.1:12; t=interp1(hours,temps,h,spline); (直接输出数据将是很多的) plot(hours,temps,+,h,t,hours,temps,r:) % 作图 xlabel(Hour),ylabel(Degrees Celsius)2、二维插值 对二维插值问题，Matlab分别给出了针对插值基点为网络节点的插值函数interp2，以及针对插值基点为散乱节点的插值函数griddata。现分别介绍如下：2.1插值基点为网格节点问题：已知个节点：且。求点处的插值。对上述问题，Matlab中提供了二维插值函数interp2，其基本格式为：cz=interp2(x,y,z,cx,cy,method)其中向量和是自变量。的分量值必须是单调递增的。是维向量，指明所给数据网格点的横坐标；是维向量，指明所给数据网格点的纵坐标；是维矩阵，标明相应于所给数据网格点的函数值。向量cx,cy是给定的网格点的横坐标和纵坐标，指明函数cz=interp2(x,y,z,cx,cy,method)返回在网格(cx,cy)处的函数值。cx与cy应是方向不同的向量。即一个是行向量；另一个是列向量。Method为可选参数，可从以下四个值中任选一个：nearest-最近邻点插值linear-线性插值spline-三次样条插值cublic-三次插值linear是缺省值例4 测的平板表面往各点处的温度分别为：作出平板表面温度分布曲面。解：1. 先在二维坐标画出原始数据，看一下该数据的粗糙程度。输入以下命令：x=1:5;y=1:3;temps=82 81 80 82 84;79 63 61 65 81;84 84 82 85 86;mesh(x,y,temps)画出粗糙的温度分布曲面图2. 在x,y方向上每隔0.2个单位的地方进行插值，以平滑数据。再输入以下命令：xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi,yi,cubic)mesh(xi,yi,zi)画出插值后的温度分布曲面图2. 2插值基点为散乱节点问题：已知个节点：，求点处的插值对上述问题，Matlab提供了插值函数griddata。其格式为:cz=griddata(x,y,z,cx,cy,method)其中x,y,z均是n维向量，指明所给数据点的横坐标纵坐标和竖坐标。向量cx,cy是给定的网格点的横坐标和纵坐标，指明函数cz=griddata(x,y,z,cx,cy,method)返回在网格（cx,cy）处的函数值。cx与cy应是方向不同的向量。即一个是行向量；另一个是列向量。Method为可选参数，可从以下四个值中任选一个：nearest-最近邻点插值linear-线性插值cublic-三次插值v4-Matlab中所提供的插值方法linear是缺省值插值函数e01sef和e01sff针对散乱基点的二维插值问题，Matlab还提供了基于修正的Shepherd插值法的函数e01sef和e01sff这两个函数必须同时使用。首先使用e01sef求出参数rnw和fnodes，再将它们用于函数e01sff中，求得插值点的坐标。函数的格式如下：fnodes,minnq,rnw,rnq,ifail=e01sef(x,y,z)cz(i,j),ifail=e01sff(x,y,z,rnw,fnodes,cx(i),cy(j)其中x,y,z均是n维向量，指明所给数据点的横坐标，纵坐标和竖坐标。向量cx,cy是给定的网格点的横坐标和纵坐标,cz是矩阵，其行数和列数分别等于cx和cy的维数。函数e01sff求出一个插值点(cx(i),cy(j)所对应的函数值cz(i,j).例5在某海域测得一些点(x,y)处的水深z由下表给出，在矩形区域(75,200)*(-50,150)内画出海底曲面的图形x129140103.588185.5195105157.5107.57781162162117.5y7.5141.52314722.5137.585.5-6.5-81356.5-66.584-33.5z48686889988949该问题的m文件hd1.m如下：%数据输入x=129 140 103.5 88 185.5 195 105 157.5 107.5 77 81 162 162 117.5 ;y=7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5;z=-4 -8 -6 -8 -6 -8 -8 -9 -9 -8 -8 -9 -4 -9;%插值cx=75:0.5:200;cy=-50:0.5:150;cz=griddata(x,y,z,cx,cy,cubic);%作表面图mesh(cx,cy,cz)三、 数据拟合曲线拟合问题是指：已知平面上个点互不相同。寻求函数，使在某种准则下与所有数据点最为接近，即曲线拟合得最好。线性最小二乘法是解决曲线拟合最常用的方法，其基本思路是，令 (1) 其中是事先选定的一组函数，是待定系数。寻求使（2）最小，称为最小二乘准则。这种拟合方法称为线性最小二乘拟合，其结果为 (2)1、系数的确定 为求使达到最小，只需利用极值的必要条件，得到关于的线性方程组线性最小二乘拟合可由解超定方程： (3) 来求得。其中,当线性无关时，可逆，方程组(3)有唯一解。在Matlab中，此解为：注：对于方程组Ax=b，A为nm矩阵，如果A列满秩，且nm。则方程组没有精确解，此时称方程组为超定方程组。线性超定方程组经常遇到的问题是数据的曲线拟合。对于超定方程，在MATLAB中，利用左除命令（x=Ab）来寻求它的最小二乘解；还可以用广义逆来求，即x=pinv(A),所得的解不一定满足Ax=b，x只是最小二乘意义上的解。左除的方法是建立在奇异值分解基础之上，由此获得的解最可靠；广义逆法是建立在对原超定方程直接进行householder变换的基础上，其算法可靠性稍逊与奇异值求解，但速度较快。2、 的选取 面对一组数据，用线性最小二乘法作曲线拟合时，关键是恰当地选取。如能通过机理分析，知道和之间应该有什么样的函数关系，则就容易确定。若无法知道和之间的关系，通常可以将数据作图，直观地判断应该用什么样的曲线去作拟合，常用的曲线有直线，多项式，双曲线，指数曲线等。实际操作中可以在直观判断的基础上，选几种曲线分别作拟合，然后比较，看哪条曲线的最小二乘指标最小。 特别地，若令 (4)则线性最小二乘拟合称为多项式拟合。3、 多项式拟合假设有已知数据。现求作一个不超过次的多项式，使得取最小，亦即 (5)取最小。不难看出，以上多项式最小二乘拟合问题就是求解关于的超定方程组 (6)的最小二乘解问题。把当作变量，上述方程组的矩阵法记为这是一个超定方程组。根据定理1，它所对应的正规方程组为记,上述方程组可改写为 (7)通过解正规方程组(7)便可解出，从而确定出拟合多项式。多项式拟合的一般方法可归纳为：(1) 根据具体问题，确定拟合多项式的次数;(2) 由公式,计算出与;(3) 写出正规方程组(7)式；(4) 解正规方程组,求出(5) 写出拟合多项式。四、 用MATLAB解曲线拟合问题在Matlab中实现最小二乘法拟合通常可以采用如下两种途径：(1) 利用polyfit函数进行多项式拟合。(2) 利用矩阵除法解决复杂型函数的拟合。 线性最小二乘拟合 在Matlab的线性最小二乘拟合中，用的较多的是多项式拟合，其命令为：其中多项式在处的值可用以下命令计算: y=polyval(A,x) 例1 设,用最小二乘法拟合如下的数据0.51.01.52.02.53.01.752.453.814.808.008.60此例用polyfit功能函数进行拟合。解：在命令窗中输入： x=0.5 1.0 1.5 2.0 2.5 3.0;(或者：x=0.5:0.5:3.0) y=1.75 2.45 3.81 4.80 8.00 8.60; a=polyfit(x,y,2) a = 0.4900 1.2501 0.8560 x1=0.5:0.05:3.0; y1=polyval(a,x1); plot(x,y,)hold onplot(x1,y1,-r)例2 对下面一组数据作二次多项拟合:00.10.20.30.40.50.60.70.80.91-0.447