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),y2,p2(7)=,x0=1,x1=4,x2=9,y0=1,y1=2,y2=3,(14)(19),(74)(79),1,+,(41)(49),(71)(79),2,+,(91)(94),(71)(74),3,=2.7,p2(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 xn,y=f(x),y=p(x),a,b,分段线性插值法,从舍入误差来看，高次插值误差的传播也较为严重，在一个节点上产生的舍入误差会在计算中不断扩大，并传播到其它节点上。因此，次数太高的高次插值多项式并不实用，因为节点数增加时，计算量增大了，但插值函数的精度并未提高。为克服在区间上进行高次插值所造成的龙格现象，采用分段插值的方法，将插值区间分成若干个小的区间，在每个小区间进行线性插值，然后相互连接，用连接相邻节点的折线逼近被插函数，这种把插值区间分段的方法就是分段线性插值法。,分段线性插值法,分段线性插值法在几何上就是用折线替代曲线,例:已知f(x)在四个节点上的函数值如下表所示,30 45 60 90,1,求f(x)在区间30,90上的分段连续线性插值函数S(x),解:将插值区间30,90分成连续的三个小区间 30,45,45,60,60,90 则S(x)在区间30,45上的线性插值为,S(x)在区间45,60上的线性插值为,S(x)在区间60,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),nearest 最邻近插值；linear 线性插值；spline 三次样条插值；cubic 立方插值；缺省时 分段线性插值,注意：所有的插值方法都要求x是单调的，并且xi不能够超过x的范围,nearest 最邻近插值；linear 线性插值；spline 三次样条插值；cubic 立方插值；缺省时 分段线性插值,nearest 最邻近插值；linear 线性插值；spline 三次样条插值；cubic 立方插值；缺省时 分段线性插值,例：从1点12点的11小时内，每隔1小时测量一次温度，测得的温度的数值依次为：5，8，9，15，25，29，31，30，22，25，27，24试估计每隔1/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),例 已知飞机下轮廓线上数据如下，求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的范围,例：测得平板表面35网格点处的温度分别为：82 81 80 82 84 79 63 61 65 81 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),nearest最邻近插值linear 双线性插值cubic 双三次插值v4-MATLAB提供的插值方法缺省时,双线性插值,要求cx取行向量，cy取为列向量,例：在某海域测得一些点(x,y)处的水深z由下表给出，船的吃水深度为5英尺，在矩形区域（75，200）（-50，150）里的哪些地方船要避免进入,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),数据拟合,对于已给一批实测数据，由于实测方法、实验环境等一些外界因素的影响，不可避免地会产生随机干扰和误差。我们自然希望根据数据分布的总趋势去剔除观察数据中的偶然误差，这就是所谓的数据修匀（或称数据平滑）问题。,换句话说:求一条曲线,使数据点均在离此曲线的上方或下方不远处,所求的曲线称为拟合曲线,它既能反映数据的总体分布,又不至于出现局部较大的波动,更能反映被逼近函数的特性,使求得的逼近函数与已知函数从总体上来说其偏差按某种方法度量达到最小,这就是曲线拟合。,数据拟合,仍然是已知 x1 xm;y1 ym,求一个简单易算的近似函数 f(x)来拟合这些数据。yi 本身是测量值，不准确，即 yi f(xi)；这时没必要取 f(xi)=yi,而要使 i=f(xi)yi 总体上尽可能地小。这种构造近似函数 的方法称为曲线拟合，f(x)称为拟合函数。拟合模型的分类 直线拟合曲线拟合,直 线 拟 合 问 题 引 例,求60C时的电阻R,设 R=at+ba,b为待定系数,曲 线 拟 合 问 题 引 例,求血药浓度随时间的变化规律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),可化为线性拟合的非线性拟合,有些非线性拟合曲线可以通过适当的变量替换转化为线性曲线，从而用线性拟合进行处理，对于一个实际的曲线拟合问题，一般先按观测值在直角坐标平面上描出散点图，看一看散点的分布同哪类曲线图形接近，然后选用相接近的曲线拟合方程。再通过适当的变量替换转化为线性拟合问题，按线性拟合解出后再还原为原变量所表示的曲线拟合方程。下表列举了几类经适当变换后化为线性拟合求解的曲线拟合方程及变换关系,曲线拟合方程 变换关系 变换后线性拟合方程,最小二乘法的优劣,我们已经讨论了最小二乘意义下的曲线拟合问题,由于方程比较简单,实际中应用广泛,特别是因为任何连续函数至少在一个较小的邻域内可以用多项式任意逼近,因此用多项式作数据拟合,有它的特殊重要性。从而在许多实际问题中,不论具体函数关系如何,都可用多项式作近似拟合,但用多项式拟合时,当n较大时(n7),其法方程的系数矩阵的条件数一般较大,所以往往是病态的,因而给求解工作带来了困难。近年来,产生一些直接解线性最小二乘问题的新方法，例如正交三角化方法。,用MATLAB作线性最小二乘拟合,1.作多项式f(x)=a1xm+amx+am+1拟合,可利用已有程序:,a=polyfit(x,y,m),2.多项式在x处的值y可用以下命令计算：y=polyval（a，x）,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：1949年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 建模分析我国人口增长的规律,预报1999年我国人口数。,模型一：假设人口随时间线性地增加,模型：参数估计观测值的模型：拟合的精度:误差平方和。可以算出：a=-283.2320 b=0.1480模型：y=1.93+0.146 x,则可看成是线性方程,用 polyfit命令计算得：,模型二：指数增长模型,可变为,Y,A,=,+,BX,a=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 人口数 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 误 差 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 误差-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：已知观测数据点如表所示,分别用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)grid,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),例2 用切削机床进行金属品加工时,为了适当地调整机床,需要测定刀具的磨损速度.在一定的时间测量刀具的厚度,得数据如表所示:,解:描出散点图,在命令窗口输入:,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.3804,hold on,plot(t,y1),hold off,a=polyfit(t,y,1),y1=-0.3012*t+29.3804,在实际应用中常见的拟合曲线,直线,多项式,双曲线(一支),指数曲线,用MATLAB作非线性最小二乘法拟合,两个求非线性最小二乘拟合的函数：lsqcurvefit、lsqnonlin。相同点和不同点：两个命令都要先建立M-文件fun.m，定义函数f(x)，但定义f(x)的方式不同,非线性曲线拟合:lsqcurvefit,功能:,x=lsqcurvefit(fun,x0,xdata,ydata),x,resnorm=lsqcurvefit(fun,x0,xdata,ydata),根据给定的数据 xdata,ydata(对应点的横,纵坐标),按函数文件 fun 给定的函数,以x0为初值作最小二乘拟合,返回函数 fun中的系数向量x和残差的平方和resnorm.,基本概念,稀疏矩阵(sparse matrix)：矩阵中非零元素的个数远远小于矩阵元素的总数），并且非零元素的分布没有规律。与之相区别的是，如果非零元素的分布存在规律（如上三角矩阵、下三角矩阵、对称矩阵），则称该矩阵为特殊矩阵。稠密矩阵:非0元素占所有元素比例较大的矩阵。若n阶矩阵A的行列式不为零，即|A|0，则称A为非奇异矩阵，否则称A为奇异矩阵。如果A为n阶对称正定矩阵。把矩阵A的行换成相应的列，得到的新矩阵称为A的转置矩阵，记作A。AA=E（E为单位矩阵，）或AA=E，则n阶实矩阵 A称为正交矩阵。,矩阵分解,矩阵分解通过将复杂矩阵表示成形式简单或具有良好数学性质（统称为简单矩阵）的组合，以便于理论分析或数值计算。通常矩阵分解将复杂矩阵分解为几个简单矩阵的乘积。求解线性方程组不可避免地要用到矩阵分解的概念。MATLAB中，线性方程组的求解主要用到三种基本的矩阵分解，即对称正定矩阵的cholesky分解、一般方程的gaussian消去法和矩阵的正交分解。这三种分解由函数chol、lu和qr完成。,矩阵分解之LU分解,矩阵的三角分解又称LU分解，它的目的是将一个矩阵分解成一个下三角矩阵L和一个上三角矩阵U的乘积，即A=LU。Matlab使用函数lu实现LU分解，其格式为：L,U=lu(A)其中U为上三角阵，L为下三角阵或其变换形式，满足LU=A。L,U,P=lu(A)U为上三角阵，L为下三角阵，P为单位矩阵的行变换矩阵，满足LU=PA。例:A=1 2 3;4 5 6;7 8 9;L,U=lu(A)L,U,P=lu(A),矩阵分解之Cholesky分解,对于正定矩阵A，则存在一个实的非奇异上三角阵R，满足R*R=A，称为Cholesky分解。Matlab使用函数chol实现Cholesky分解，其格式为：R=chol(A)若A非正定，则产生错误信息。R,p=chol(A)不产生任何错误信息，若A为正定阵，则p=0，R与上相同；若A非正定，则p为正整数，R是有序的上三角阵。例：A=pascal(4);R,p=chol(A),矩阵分解之QR分解,将矩阵A分解成一个正交矩阵Q与一个上三角矩阵R的乘积A=QR，称为QR分解。Matlab使用函数qr实现QR分解，其格式为：Q,R=qr(A)Q,R,E=qr(A)求得正交矩阵Q和上三角阵R，E为单位矩阵的变换形式，R的对角线元素按大小降序排列，满足AE=QR。Q,R=qr(A,0)例：A=1 2 3;4 5 6;7 8 9;10 11 12;Q,R=qr(A),线性方程的组的求解,我们将线性方程的求解分为两类：一类是方程组求唯一解或求特解，另一类是方程组求无穷解即通解。可以通过系数矩阵的秩来判断：若系数矩阵的秩r=n(n为方程组中未知变量的个数)，则有唯一解；若系数矩阵的秩rn，则可能有无穷解；求秩函数：RANK(A),线性方程的组的求解,求线性方程组的唯一解或特解（直接法）利用矩阵除法求线性方程组AX=b的特解利用LU、QR和cholesky分解求方程组的解求线性方程组的通解线性方程组的无穷解=对应齐次方程组的通解+非齐次方程组的一个特解,求线性方程组的唯一解或特解,利用矩阵除法求线性方程组AX=b的特解当系数矩阵A为N*N的方阵时，MATLAB会自行用高斯消去法求解线性代数方程组。若右端项b为N*1的列向量，则x=Ab可获得方程组的数值解x（N*1的列向量）若右端项b为N*M的矩阵，则x=Ab可同时获得同一系数矩阵A、M个方程组数值解x（为N*M的矩阵），即x(:,j)=Ab(:,j)，j=1,2,M。,求线性方程组的唯一解或特解,用函数rref求解C=A,B%由系数矩阵和常数列构成增广矩阵CR=rref(C)%将C化成行最简行则R的最后一列元素就是所求之解。,例1：求下列方程组的解。,A=5 6 0 0 0 1 5 6 0 0 0 1 5 6 0 0 0 1 5 6 0 0 0 1 5;B=1 0 0 0 1;R_A=rank(A)%求秩X=AB%求解,例2：求下列方程组的解。,A=1 1-3-1;3-1-3 4;1 5-9-8;B=1 4 0;R_A=rank(A)%求秩X=AB%由于系数矩阵不满秩,该解法可能存在误差X=0 0-0.5333 0.6000%一个特解近似值,A=1 1-3-1;3-1-3 4;1 5-9-8;B=1 4 0;C=A,B;%构成增广矩阵R=rref(C)若用rref求解，则比较精确：,LU分解求方程组的解,LU分解又称Gauss消去分解，可把任意方阵分解为下三角矩阵和上三角矩阵的乘积。即A=LU，L为下三角阵，U为上三角阵。则：A*X=b L*U*X=b所以X=U(Lb)这样可以大大提高运算速度。若采用l，u，p=lu(A)，则Ax=b的解为 x=u(l(p*b),利用LU分解求方程组的解例3:求下列方程组的一个特解。解：A=4 2-1;3-1 2;11 3 0;B=2 10 8;L,U=lu(A)X=U(LB),Cholesky分解求方程组的解,若A为对称正定矩阵，则Cholesky分解可将矩阵A分解成上三角矩阵和其转置的乘积，即：A=R*R 其中R为上三角阵。方程A*X=b变成R*R*X=b所以X=R(Rb),QR分解求方程组的解,对于任何长方矩阵A，都可以进行QR分解，其中Q为正交矩阵，R为上三角矩阵的初等变换形式，即：A=QR方程A*X=b变形成QRX=b所以X=R(Qb)若采用Q，R，E=qr(X)，则Ax=b的解为 x=E(R(Qb),求线性齐次方程组的通解,在Matlab中，函数null用来求解零空间，即满足AX=0的解空间，实际上是求出解空间的一组基(基础解系)。其格式为：z=null(A)z的列向量为方程组的正交规范基z=null(A,r)z的列向量是方程AX=0的有理基,例4:求下列方程组的通解。A=1 2 2 1;2 1-2-2;1-1-4-3;format rat%指定有理式格式输出B=null(A,r)%求解空间的有理基,对于线性方程组Ax=b情形1：m=n（恰定方程组）情形2：mn（超定方程），多用于曲线拟合。注：m为表达式个数，n为未知数个数,求非齐次线性方程组的通解,非齐次线性方程组需要先判断方程组是否有解，若有解，再去求通解。因此，步骤为：第