MATLAB的数值分析.ppt

第十章MATLAB的数值分析制作:陈学明数据拟合在第53页幻灯片集中趋势的测定在统计研究中，需要搜集大量数据并对其进行加工整理，对这些数据进行整理之后发现：大多数情况下数据都会呈现出一种钟形分布，即各个变量值与中间位置的距离越近，出现的次数越多；与中间位置距离越远，出现的次数越少，从而形成了一种以中间值为中心的集中趋势。这个集中趋势是现象共性的特征，是现象规律性的数量表现。集中趋势的描述一均值函数算术平均数 mean(x)调和平均数harmmean(x)几何平均数 geomean(x)切尾均值trimmean(x，percent)集中趋势的描述二中位数（中位次数）函数中位数是指全体数值按大小排列后，位于中间的数值。如果参数集合中包含有偶数个数字，中位数函数将返回位于中间的两个数的平均值。Median(x)集中趋势的描述三众数函数众数是一组数列中出现次数最多的数值，众数函数返回某一数组或数据区域中出现频率最多的数值。最大（小）值函数最大（小）值函数可以返回数据集中的最大（小）数值。range(X)Max(x)-min(x)三种平均数的特点众数是一组数据中出现次数最多的变量值，它用于对分类数据的概括性度量，其特点是不受极端值的影响，但它没有利用全部数据信息，而且还具有不惟一性。一组数据可能有众数，也可能没有众数；可能有一个众数，也可能有多个众数。中位数是一组数据按大小顺序排序后处于中间位置上的变量，它主要用于对顺序数据的概括性度量。均值是一组数据的算术平均，它利用了全部数据信息，是概括一组数据最常用的一个值。离中趋势的测定在研究现象总体标志的一般水平时，不仅要研究总体标志的集中趋势，还要研究总体标志的离中趋势，如研究价格背离价值的平均程度。研究离中趋势可以通过计算标志变异指标来进行。标志变异指标是同统计平均数相联系的一种综合指标，用于度量随机变量在取值区间内的分布情况，主要有平均差、标准差、方差、四分位数、百分位数等。样本标准差样本标准差函数用来估算样本的标准偏差，反映相对于平均值（mean）的离散程度，计算样本标准差采用不偏估计式（亦即自由度n-1），其计算公式为:std(X)总体标准差总体标准差函数返回以参数形式给出的整个样本总体的标准偏差，反映相对于平均值（mean）的离散程度。计算总体标准差使用整个总体的变量，通常采用偏性估计式（亦即自由度为n），其计算公式为四分位数四分位数是将中值的前后两部分数值再等分为二，以数值小的一端算起，前半部的分区点称为第1四分位数，后半部的分区点称为第3四分位数，而中值即为第2四分位数。四分位数通常用于在销售额和测量值数据集中对总体进行分组。w为数据序列，Q1为上四分位值，Q3为下四分位值，计算如下：Q1=prctile(w,25);Q3=prctile(w,75);prctile()函数实现计算样本的百分位数功能分布形态的测定只用集中趋势和离中趋势来表示所有数据，难免不够准确。分析总体次数的分布形态有助于识别整个总体的数量特征。总体的分布形态可以从两个角度考虑，一是分布的对称程度，另一个是分布的高低。前者的测定参数称为偏度或偏斜度，后者的测定参数称为峰度。峰度是掌握分布形态的另一指标，它能描述分布的平缓或陡峭程度。如果峰度数值等于零，说明分布为正态；若峰度数值大于零，说明分布呈陡峭状态；若峰度数值小于零，说明分布形态趋于平缓。偏度函数偏度函数返回分布的偏斜度。偏斜度反映以平均值为中心的分布的不对称程度。正偏斜度表示不对称边的分布更趋向正值（偏右），负偏斜度表示不对称边的分布更趋向负值（偏左）。其计算公式为skewness(X)峰度函数峰度函数返回数据集的峰值，表示次数分布高峰的起伏状态。峰值反映与正态分布相比某一分布的尖锐度或平坦度。正峰值表示相对尖锐的分布，负峰值表示相对平坦的分布。其计算公式为：kurtosis(X)插值与拟合插值与拟合属数值分析中函数逼近内容。在数学建模中，插值与拟合是一种常用的数据分析手段，被公认为建模中的十大算法之一。插值问题与拟合问题矿井中某处的瓦斯浓度y与该处距地面的距离x有关，现用仪器测得从地面到井下500米每隔50米的瓦斯浓度数据(xi,yi)(i=0,1,10)，根据这些数据完成下列工作:(1)寻找一个函数，要求由此函数可近似求得从地面到井下500米之间任意点处的瓦斯浓度；(2)估计井下600米处的瓦斯浓度。第一个问题可归结为“已知函数在x0,x1,xn处的值，求函数在区间x0,xn内其它点处的值”，这种问题适宜用插值方法解决。插值问题可描述为：已知函数在x0,x1,xn处的值y0,y1,yn，求函数p(x)，使p(xi)=yi。但对第二个问题不宜用插值方法,因为600米已超出所给数据范围，用插值函数外推插值区间外的数据会产生较大的误差。第一个问题可归结为第一个问题可归结为“已知函数在已知函数在x0,x1,xn处的值，求函数在区间处的值，求函数在区间x0,xn内其它点处的值内其它点处的值”，这种问题适宜用插值方法解决。这种问题适宜用插值方法解决。插值问题可描述为：已知函数在插值问题可描述为：已知函数在x0,x1,xn处的值处的值y0,y1,yn，求函数，求函数p(x)，使，使p(xi)=yi。但对第二个问题不宜用插值方法但对第二个问题不宜用插值方法,因为因为600米已超出米已超出所给数据范围，用插值函数外推插值区间外的数据所给数据范围，用插值函数外推插值区间外的数据会产生较大的误差。会产生较大的误差。解决第二个问题的常用方法是，根据地面到井下解决第二个问题的常用方法是，根据地面到井下 500 处的处的数据求出瓦斯浓度与地面到井下距离数据求出瓦斯浓度与地面到井下距离x之间的近似函数关之间的近似函数关系系f(x),由由f(x)求井下求井下600米处的瓦斯浓度。米处的瓦斯浓度。插值函数过已知点，拟合函数不一定过已知点。通插值函数过已知点，拟合函数不一定过已知点。通常常,插值主要用于求函数值，而拟合的主要目的是求插值主要用于求函数值，而拟合的主要目的是求函数关系。当然，某些问题既可以用插值也可以用函数关系。当然，某些问题既可以用插值也可以用拟合。拟合。插值方法-概述为什么需要插值？(1)函数关系y=f(x)没有明确的表达式(2)y=f(x)表达式复杂，不便于研究和使用 一般地称g(x)为逼近函数，f(x)为被逼近函数。如果要求构造的函数g(x)取给定的离散数据，则称g(x)为f(x)的插值函数。当g(x)为代数多项式时，与之相应的插值逼近称为代数多项式插值。插值多项式存在且唯一插值多项式存在且唯一一维插值若已知的数据集是平面上的一组离散点集(x,y)，则其相应的插值就是一维插值。常用插值方法：拉格朗日插值（lagrange插值）分段线性插值三次样条插值拉格朗日插值在数值分析中，拉格朗日插值法是以法国十八世纪数学家约瑟夫路易斯拉格朗日命名的一种多项式插值方法。许多实际问题中都用函数来表示某种内在联系或规律，而不少函数都只能通过实验和观测来了解。如对实践中的某个物理量进行观测，在若干个不同的地方得到相应的观测值，拉格朗日插值法可以找到一个多项式，其恰好在各个观测的点取到观测到的值。这样的多项式称为拉格朗日（插值）多项式。数学上来说，拉格朗日插值法可以给出一个恰好穿过二维平面上若干个已知点的多项式函数。拉格朗日插值法最早被英国数学家爱德华华林于1779年发现1，不久后（1783年）由莱昂哈德欧拉再次发现。1795年，拉格朗日在其著作师范学校数学基础教程中发表了这个插值方法，从此他的名字就和这个方法联系在一起拉格朗日(Lagrange)插值已知函数f(x)在n+1个点x0,x1,xn处的函数值为y0,y1,yn。求一n次多项式函数Pn(x)，使其满足：Pn(xi)=yi,i=0,1,n.解决此问题的拉格朗日插值多项式公式如下其中Li(x)为n次多项式：拉格朗日(Lagrange)插值特别地:两点一次(线性)插值多项式三点二次(抛物)插值多项式线性插值线性插值的几何意义线性插值的几何意义:用通过用通过A、B两点的直线近似地代替两点的直线近似地代替曲线曲线 y=f(x)。由解析几何知道由解析几何知道,这条直线用这条直线用点斜式表示为点斜式表示为例：例：已知已知y=f(x)的函数表的函数表 求线性插值多项式求线性插值多项式,并计算并计算x=1.5 的值的值X 1 3 y 1 2解解:由线性插值多项式公式得由线性插值多项式公式得抛物插值抛物插值又称二次插值，它也是常用的代数插值之一。设已知f(x)在三个互异点x0，x1，x2的函数值y0，y1，y2，要构造次数不超过二次的多项式使满足二次插值条件：使满足二次插值条件：例：已知x=1,4,9的平方根值,用抛物插值公式,求(x0 x1)(x0 x2)(xx1)(xx2)y0+(x1x0)(x1x2)(xx0)(xx2)y1+(x2x0)(x2x1)(xx0)(xx1)y2p2(7)=x0=1,x1=4,x2=9y0=1,y1=2,y2=3(14)(19)(74)(79)1+(41)(49)(71)(79)2+(91)(94)(71)(74)3=2.7p2(x)=2.6458拉格朗日插值多项式两个插值点可求出一次插值多项式两个插值点可求出一次插值多项式,而三个插值点可而三个插值点可求出二次插值多项式。插值点增加到求出二次插值多项式。插值点增加到n+1个时个时,也就是也就是通过通过n+1个不同的已知点来构造一个次数为个不同的已知点来构造一个次数为n的代数的代数多项式多项式P(x)。拉格朗日插值多项式结构对称，使用方便。但由于是拉格朗日插值多项式结构对称，使用方便。但由于是用基函数构成的插值，这样要增加一个节点时，所有用基函数构成的插值，这样要增加一个节点时，所有的基函数必须全部重新计算，不具备承袭性，还造成的基函数必须全部重新计算，不具备承袭性，还造成计算量的浪费。计算量的浪费。多数情况下，Lagrange插值法效果是不错的，但随着节点数n的增大，Lagrange多项式的次数也会升高，可能造成插值函数的收敛性和稳定性变差。如龙格（Runge）现象。高次插值中的Runge现象通常选用多项式作为插值函数。在研究插值问题的初期，所有人都认为插值多项式的次数越高，插值精度越高。Runge通过对一个例子的研究发现，上述结论仅仅在插值多项式的次数不超过七时成立；插值多项式的次数超过七时，插值多项式会出现严重的振荡现象，称之为Runge现象。避免Runge现象的常用方法是：将插值区间分成若干小区间，在小区间内用低次(二次，三次)插值，即分段低次插值，如样条函数插值。考察函数考察函数 右图给出了右图给出了和和 的图像的图像,当当n增大时增大时,在两端在两端会发出激烈的振荡会发出激烈的振荡,这就是所谓龙格现这就是所谓龙格现象。该现象表明象。该现象表明,在在大范围内使用高次大范围内使用高次插值插值,逼近的效果往逼近的效果往往是不理想的往是不理想的 高次多项式插值产生的不收敛现象高次多项式插值产生的不收敛现象-Runge 现象现象插值多项式的误差在插值区间a,b上用插值多项式p(x)近似代替f(x),除了在插值节点xi上没有误差外，在其它点上一般是存在误差的。若记若记 R(x)=f(x)-p(x)则则 R(x)就是用就是用 p(x)近似代替近似代替 f(x)时时的截断误差的截断误差,或称插值余项。或称插值余项。x0 x1 xixi+1 xn-1 xny=f(x)y=p(x)ab分段线性插值法从舍入误差来看，高次插值误差的传播也较为严重，在一个节点上产生的舍入误差会在计算中不断扩大，并传播到其它节点上。因此，次数太高的高次插值多项式并不实用，因为节点数增加时，计算量增大了，但插值函数的精度并未提高。为克服在区间上进行高次插值所造成的龙格现象，采用分段插值的方法，将插值区间分成若干个小的区间，在每个小区间进行线性插值，然后相互连接，用连接相邻节点的折线逼近被插函数，这种把插值区间分段的方法就是分段线性插值法。分段线性插值法分段线性插值法在几何上就是用折线替代曲线例例:已知已知f(x)f(x)在四个节点上的函数值如下表所示在四个节点上的函数值如下表所示 30 45 60 901求求f(x)f(x)在区间在区间 30,9030,90 上的分段连续线性插值函数上的分段连续线性插值函数S(x)S(x)解解:将插值区间将插值区间 30,9030,90 分成连续的三个小区间分成连续的三个小区间 30,4530,45,45,6045,60,60,9060,90 则则S(x)在区间在区间 30,4530,45 上的线性插值为上的线性插值为 S(x)在区间在区间 45,6045,60 上的上的线性插值为线性插值为 S(x)S(x)在区间在区间 60,9060,90 上的线性插值上的线性插值为为 将各小区间的将各小区间的线性插值函数连接在一起，得线性插值函数连接在一起，得 样条函数插值法给定n+1个节点上的函数值可以作n次插值多项式,但当n较大时,高次插值不仅计算复杂,而且可能出现Runge现象,采用分段插值虽然计算简单、也有一致收敛性,但不能保证整条曲线在连接点处的光滑性如分段线性插值,其图形是锯齿形的折线,虽然连续,但处都是“尖点”,因而一阶导数都不存在,这在实用上,往往不能满足某些工程技术的高精度要求。如在船体、飞机等外形曲线的设计中,不仅要求曲线连续,而且要有二阶光滑度,即有连续的二阶导数这就要求分段插值函数在整个区间上具有连续的二阶导数。因此有必要寻求一种新的插值方法这就是样条函数插值法三次样条函数设函数定义在区间a,b上，给定n+1个节点和一组与之对应的函数值，若函数满足:在每个节点上满足S(xi)=f(xi)(i=0,1,n)在a,b上有连续的二阶导数在每个小区间xi,xi+1(i=0,1,n-1)上是一个三次多项式。则称S(x)为三次样条插值函数。三次样条插值函数S(x)是一个分段三次多项式。另一方面,要求分段三次多项式S(x)及其导数S(x)和S(x)在整个插值区间a,b上连续,则要求它们在各个子区间的连接点上连续。三次样条函数用三次样条绘制的曲线不仅有很好的光滑度，而且当节点逐渐加密时，其函数值在整体上能很好地逼近被插函数，相应的导数值也收敛于被插函数的导数，不会发生龙格现象。因此三次样条在计算机辅助设计中有广泛的应用。几种插值方法的比较插值法中的拉格朗日插值多项式是研究数值微积分与微分方程数值解的重要工具。牛顿插值多项式是拉格朗日插值多项式的变形，具有承袭性，比拉格朗日插值多项式节省计算量。分段低次多项式插值由于具有良好的稳定性与收敛性，且算法简单，便于应用。特别是应用广泛的三次样条插值，不但有较好的稳定性和收敛性，而且具有较好的光滑性，从而满足了许多实际问题的要求。一维插值函数yi=interp1(x，y，xi，method)xi处的插处的插值结果值结果插值节点插值节点被插值点被插值点插值方法插值方法nearest 最邻近插值；最邻近插值；linear 线性插值；线性插值；spline 三次样条插值；三次样条插值；cubic 立方插值；立方插值；缺省时缺省时 分段线性插值分段线性插值 注意：所有的插值注意：所有的插值方法都要求方法都要求x是单调的，是单调的，并且并且xi不能够超过不能够超过x的的范围范围nearest 最邻近插值；最邻近插值；linear 线性插值；线性插值；spline 三次样条插值；三次样条插值；cubic 立方插值；立方插值；缺省时缺省时 分段线性插值分段线性插值nearest 最邻近插值；最邻近插值；linear 线性插值；线性插值；spline 三次样条插值；三次样条插值；cubic 立方插值；立方插值；缺省时缺省时 分段线性插值分段线性插值例：从例：从1 1点点1212点点的的1111小时内，每隔小时内，每隔1 1小时测量一次温小时测量一次温度，测得的温度的数值依次为：度，测得的温度的数值依次为：5 5，8 8，9 9，1515，2525，2929，3131，3030，2222，2525，2727，2424试估计每隔试估计每隔1/101/10小小时的温度值时的温度值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)xy机翼下机翼下轮廓线轮廓线例例 已知飞机下轮廓线上数据如下，求已知飞机下轮廓线上数据如下，求x每改变每改变0.1时的时的y值值x0=0 3 5 7 9 11 12 13 14 15;y0=0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6;x=0:0.1:15;%插值及作图在下页插值及作图在下页plot(x0,y0,kp,x,y1,r)plot(x0,y0,kp,x,y2,r)plot(x0,y0,kp,x,y3,r)y1=interp1(x0,y0,x,nearest);y2=interp1(x0,y0,x);y3=interp1(x0,y0,x,spline);二维插值函数（网格点）z=interp2(x0,y0,z0,x,y,method)被插值点的函数值被插值点的函数值插值节点插值方法nearest 最邻近插值；最邻近插值；linear 双线性插值；双线性插值；cubic 双三次插值；双三次插值；缺省时缺省时 双线性插值双线性插值.要求要求x0,y0单调；单调；x，y可取可取为矩阵，或为矩阵，或x取行取行向量，向量，y取为列向量，取为列向量，x,y的值分别不能超出的值分别不能超出x0,y0的范围的范围例：测得平板表面例：测得平板表面3 35 5网格点处的温度分别为：网格点处的温度分别为：82 81 80 82 84 82 81 80 82 84 79 63 61 65 81 79 63 61 65 81 84 84 82 85 86 84 84 82 85 86 试作出平板表面的温度分布曲面试作出平板表面的温度分布曲面z=f(x,y)的图形的图形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)1.先在三维坐标画出原始数据，画出粗糙的温度分布曲线图先在三维坐标画出原始数据，画出粗糙的温度分布曲线图.xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi,yi,cubic);mesh(xi,yi,zi)画出插值后的温度分布曲面图画出插值后的温度分布曲面图.2以平滑数据以平滑数据,在在 x、y方向上每隔方向上每隔0.2个单位的地方进行插值个单位的地方进行插值.通过此例对最近邻点插值、双线性插值方法和双三次插值方法的插通过此例对最近邻点插值、双线性插值方法和双三次插值方法的插值效果进行比较值效果进行比较x=0:400:5600;y=0:400:4800;z=370 470 550 600 670 690 670 620 580 450 400 300 100 150 250;.510 620 730 800 850 870 850 780 720 650 500 200 300 350 320;.650 760 880 970 1020 1050 1020 830 900 700 300 500 550 480 350;.740 880 1080 1130 1250 1280 1230 1040 900 500 700 780 750 650 550;.830 980 1180 1320 1450 1420 1400 1300 700 900 850 840 380 780 750;.880 1060 1230 1390 1500 1500 1400 900 1100 1060 950 870 900 930 950;.910 1090 1270 1500 1200 1100 1350 1450 1200 1150 1010 880 1000 1050 1100;.950 1190 1370 1500 1200 1100 1550 1600 1550 1380 1070 900 1050 1150 1200;.1430 1430 1460 1500 1550 1600 1550 1600 1600 1600 1550 1500 1500 1550 1550;.1420 1430 1450 1480 1500 1550 1510 1430 1300 1200 980 850 750 550 500;.1380 1410 1430 1450 1470 1320 1280 1200 1080 940 780 620 460 370 350;.1370 1390 1410 1430 1440 1140 1110 1050 950 820 690 540 380 300 210;.1350 1370 1390 1400 1410 960 940 880 800 690 570 430 290 210 150;figure(1);meshz(x,y,z)xi=0:50:5600;yi=0:50:4800;figure(2)z1i=interp2(x,y,z,xi,yi,nearest);surfc(xi,yi,z1i)figure(3)z2i=interp2(x,y,z,xi,yi,linear);surfc(xi,yi,z2i)figure(4)z3i=interp2(x,y,z,xi,yi,cubic);surfc(xi,yi,z3i)二维插值函数（离散点）z=interp2(x0,y0,z0,cx,cy,method)被插值点的函数值被插值点的函数值插值节点插值方法nearestnearest最邻近插值最邻近插值linear linear 双线性插值双线性插值cubic cubic 双三次插值双三次插值v4-MATLABv4-MATLAB提供的插值方法提供的插值方法缺省时缺省时,双线性插值双线性插值要求要求cx取行向量，取行向量，cy取为列向量取为列向量 例：在某海域测得一些点例：在某海域测得一些点(x,y)处的水深处的水深z由下表给出，由下表给出，船的吃水深度为船的吃水深度为5 5英尺，在矩形区域（英尺，在矩形区域（7575，200200）（-50-50，150150）里的哪些地方船要避免进入）里的哪些地方船要避免进入1.输入插值基点数据输入插值基点数据 2.在矩形区域在矩形区域(75,200)(-50,150)进行插值。进行插值。3.作海底曲面图作海底曲面图 4.作出水深小于作出水深小于5的海域范围的海域范围,即即z=5的等高线的等高线.%程序一：插值并作海底曲面图程序一：插值并作海底曲面图x =129.0 140.0 103.5 88.0 185.5 195.0 105.5 157.5 107.5 77.0 81.0 162.0 162.0 117.5;y=7.5 141.5 23.0 147.0 22.5 137.5 85.5 -6.5 -81 3.0 56.5 -66.5 84.0 -33.5;z=4 8 6 8 6 8 8 9 9 8 8 9 4 9;x1=75:1:200;y1=-50:1:150;x1,y1=meshgrid(x1,y1);z1=griddata(x,y,z,x1,y1,v4);meshc(x1,y1,z1)%程序二：插值并作出水深小于程序二：插值并作出水深小于5的海域范围。的海域范围。x1=75:1:200;y1=-50:1:150;x1,y1=meshgrid(x1,y1);z1=griddata(x,y,z,x1,y1,v4);%插值插值z1(z1=5)=nan;%将水深大于将水深大于5的置为的置为nan，这样绘图就不会显示出来，这样绘图就不会显示出来meshc(x1,y1,z1)数据拟合对于已给一批实测数据，由于实测方法、实验环境等一些外界因素的影响，不可避免地会产生随机干扰和误差。我们自然希望根据数据分布的总趋势去剔除观察数据中的偶然误差，这就是所谓的数据修匀（或称数据平滑）问题。换句话说:求一条曲线,使数据点均在离此曲线的上方或下方不远处,所求的曲线称为拟合曲线,它既能反映数据的总体分布,又不至于出现局部较大的波动,更能反映被逼近函数的特性,使求得的逼近函数与已知函数从总体上来说其偏差按某种方法度量达到最小,这就是曲线拟合。数据拟合仍然是已知x1xm;y1ym,求一个简单易算的近似函数f(x)来拟合这些数据。yi 本身是测量值，不准确，即yif(xi)；这时没必要取f(xi)=yi,而要使i=f(xi)yi 总体上尽可能地小。这种构造近似函数的方法称为曲线拟合，f(x)称为拟合函数。拟合模型的分类直线拟合曲线拟合直线拟合问题引例温度温度t(C)20.5 32.7 51.0 73.0 95.7电阻电阻R()765 826 873 942 1032已知热敏电阻数据：已知热敏电阻数据：求求6060C C时的电阻时的电阻R R设设R=at+ba,b为待定系数为待定系数曲线拟合问题引例t(h)0.25 0.5 1 1.5 2 3 4 6 8c(g/ml)19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01已知一室模型快速静脉注射下的血药浓度数据已知一室模型快速静脉注射下的血药浓度数据(t=0注射注射300mg)求血药浓度随时间的变化规律求血药浓度随时间的变化规律c(t).曲线拟合问题最常用的解法线性最小二乘法的基本思路第一步:先选定一组函数先选定一组函数r1(x),r2(x),rm(x),mn,令令f(x)=a1r1(x)+a2r2(x)+amrm(x)（1）其中其中a1,a2,am为待定系数为待定系数 第二步:确定确定a1,a2,am的准则（最小二乘准则）：的准则（最小二乘准则）：使使n个点个点（xi,yi)与与曲线曲线 y=f(x)的距离的距离 i 的平方和最小的平方和最小 记记问题归结为，求问题归结为，求a1,a2,am 使使 J(a1,a2,am)最小最小线性最小二乘拟合函数的选取通过机理分析建立数学模型来确定f(x)将数据(xi,yi)i=1,n作图，通过直观判断确定f(x)+f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1+a2/xf=aebxf=ae-bx可化为线性拟合的非线性拟合有些非线性拟合曲线可以通过适当的变量替换转化为线性曲线，从而用线性拟合进行处理，对于一个实际的曲线拟合问题，一般先按观测值在直角坐标平面上描出散点图，看一看散点的分布同哪类曲线图形接近，然后选用相接近的曲线拟合方程。再通过适当的变量替换转化为线性拟合问题，按线性拟合解出后再还原为原变量所表示的曲线拟合方程。下表列举了几类经适当变换后化为线性拟合求解的曲线拟合方程及变换关系 曲线拟合方程曲线拟合方程 变换关系变换关系 变换后线性拟合方程变换后线性拟合方程最小二乘法的优劣我们已经讨论了最小二乘意义下的曲线拟合问题,由于方程比较简单,实际中应用广泛,特别是因为任何连续函数至少在一个较小的邻域内可以用多项式任意逼近,因此用多项式作数据拟合,有它的特殊重要性。从而在许多实际问题中,不论具体函数关系如何,都可用多项式作近似拟合,但用多项式拟合时,当n较大时(n7),其法方程的系数矩阵的条件数一般较大,所以往往是病态的,因而给求解工作带来了困难。近年来,产生一些直接解线性最小二乘问题的新方法，例如正交三角化方法。用MATLAB作线性最小二乘拟合1.作多项式作多项式f(x)=a1xm+amx+am+1拟合拟合,可可利用已有程序利用已有程序:a=polyfit(x,y,m)输入同长度输入同长度的数组的数组x，y拟合多项拟合多项式次数式次数2.2.多项式在多项式在x处的值处的值y可用以下命令计算：可用以下命令计算：y=polyval（a，x）输出拟合多项式系数输出拟合多项式系数a=a1,am,am+1（数组）数组）1）输入以下命令：）输入以下命令：x=0:0.1:1;y=-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2;A=polyfit(x,y,2)z=polyval(A,x);plot(x,y,k+,x,z,r)%作作出出数数据据点点和和拟拟合合曲曲线的图形线的图形2）计算结果：）计算结果：=-9.8108 20.1293 -0.0317用多项式拟合的命令用多项式拟合的命令例例1 1：19491949年年19941994年我国人口数据资料如下：年我国人口数据资料如下：年年 份份xi 1949 1954 1959 1964 1969 1974 1979 xi 1949 1954 1959 1964 1969 1974 1979 1984 1989 1994 1984 1989 1994 人口数人口数yi 5.4 6.0 6.7 7.0 8.1 9.1 9.8 yi 5.4 6.0 6.7 7.0 8.1 9.1 9.8 10.3 11.3 11.8 10.3 11.3 11.8 建模分析我国人口增长的规律建模分析我国人口增长的规律,预报预报19991999年我国人口年我国人口数。数。模型一：假设人口随时间线性地增加模型一：假设人口随时间线性地增加 模型：模型：参数估计观测值的模型：参数估计观测值的模型：拟合的精度拟合的精度:误差平方和。误差平方和。可以算出：可以算出：a=-283.2320 b=0.1480模型：模型：y=1.93+0.146 x 则可看成是线性方程则可看成是线性方程,用用 polyfit命令计算得：命令计算得：模型二：指数增长模型模型二：指数增长模型 可变为可变为YA=+BXa=2.33,b=0.0179则所求模型为则所求模型为:程序如下：程序如下：x=1949 1954 1959 1964 1969 1974 1979 1984 1989 1994;y=5.4 6.0 6.7 7.0 8.1 9.1 9.8 10.3 11.3 11.8;a=polyfit(x,y,1);x1=1949:5:1994;y1=a(2)+a(1)*x1;b=polyfit(x,log(y),1);y2=exp(b(2)*exp(b(1)*x1);plot(x,y,*)hold on plot(x1,y1,-r)hold on plot(x1,y2,-k)legend(原曲线原曲线,模型一曲线模型一曲线,模型二曲线模型二曲线)结论的比较如下表：结论的比较如下表：年年 份份 xi 1949 1954 1959 1964 1969 1974 1979 1984 1989 1994 xi 1949 1954 1959 1964 1969 1974 1979 1984 1989 1994 人口数人口数 yi 5.4 6.0 6.7 7.0 8.1 9.1 9.8 10.3 11.3 11.8 yi 5.4 6.0 6.7 7.0 8.1 9.1 9.8 10.3 11.3 11.8 模型一值模型一值 5.24 5.97 6.70 7.43 8.16 8.90 9.62 10.36 11.09 11.82 5.24 5.97 6.70 7.43 8.16 8.90 9.62 10.36 11.09 11.82 误误 差差 0.16 0.03 0.00 -0.43 -0.06 0.20 0.18 -0.06 0.01 -0.02 0.16 0.03 0.00 -0.43 -0.06 0.20 0.18 -0.06 0.01 -0.02 模型二值模型二值 5.55 6.06 6.62 7.23 7.90 8.64 9.44 10.31 11.26 12.31 5.55 6.06 6.62 7.23 7.90 8.64 9.44 10.31 11.26 12.31 误差误差 -0.15 -0.06 0.08 -0.23 0.20 0.46 0.36 -0.01 -0.13 -0.51 -0.15 -0.06 0.08 -0.23 0.20 0.46 0.36 -0.01 -0.13 -0.51结果分析：结果分析：（1）Q1 =0.2915 0.7437=Q2.线性模型线性模型更适合中国人口的增长。更适合中国人口的增长。（2）预报：预报：1999 年年 12.55 亿，亿，13.43 亿亿（3）统计年鉴：统计年鉴：2005 年年 13.3 亿，亿，2010 年年 14 亿亿 模型模型 I 2005 年年13.43 亿，亿，2010 年年14.16 亿亿 模型模型 II 2005 年年14.94 亿，亿，2010 年年 16.33 亿亿例例2：已知观测数据点如表所示已知观测数据点如表所示xy0-0.4470.11.9780.23.280.36.160.47.080.57.340.67.660.79.560.89.480.99.3111.2分别用分别用3次和次和6次多项式曲线拟合这些数据点次多项式曲线拟合这些数据点.x=0:0.1:1y=-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.3,11.2plot(x,y,k.,markersize,25)axis(0 1.3-2 16)p3=polyfit(x,y,3)p6=polyfit(x,y,6)编写编写Matlab程序如下程序如下:t=0:0.1:1.2s=polyval(p3,t)s1=polyval(p6,t)hold onplot(t,s,r-,linewidth,2)plot(t,s1,b-,linewidth,2)gridx=0:0.1:1y=-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.3,11.2plot(x,y,k.,markersize,25)axis(0 1.3-2 16)p3=polyfit(x,y,3)p6=polyfit(x,y,6)例例2 用切削机床进行金属品加工时用切削机床进行金属品加工时,为了适当地调整为了适当地调整机床机床,需要测定刀具的磨损速度需要测定刀具的磨损速度.在一定的时间测量在一定的时间测量刀具的厚度刀具的厚度,得数据如表所示得数据如表所示:切削时间切削时间 t/h030.0129.1228.4328.1428.0527.7627.5727.2827.0刀具厚度刀具厚度 y/cm切削时间切削时间 t/h926.81026.51126.31226.11325.71425.31524.81624.0刀具厚度刀具厚度 y/cm解解:描出散点图描出散点图,在命令窗口输入在命令窗口输入:t=0:1:16y=30.0 29.1 28.4 28.1 28.0 27.7 27.5 27.2 27.0 26.8 26.5 26.3 26.1 25.7 25.3 24.8 24.0plot(t,y,*)解解:描出散点图描出散点图,在命令窗口输入在命令窗口输入:t=0:1:16y=30.0 29.1 28.4 28.1 28.0 27.7 27.5 27.2 27.0 26.8 26.5 26.3 26.1 25.7 25.3 24.8 24.0plot(t,y,*)a=-0.3012 29.3804hold onplot(t,y1),hold offa=polyfit(t,y,1)y1=-0.3012*t+29.3804在实际应用中常见的拟合曲线直线直线多项式多项式双曲线双曲线(一支一支)指数曲线指数曲线用MATLAB作非线性最小二乘法拟合两个求非线性最小二乘拟合的函数：lsqcurvefit、lsqnonlin。相同点和不同点：两个命令都要先建立M-文件fun.m，