数学建模matlab插值与拟合ppt课件.ppt

2012数学建模培训,培训内容1.插值与拟合；2.灰色系统；3.层次分析法；(重要)4.模糊数学方法；(重要)5.时间序列与马氏链；6.多元分析与Spss；(重要)7.数学规划与Lingo；(重要)8.图论模型及其Matlab程序；,9.建模论文写作；10.建模案例讲评。培训内容基本上覆盖了数学建模中大部分常用的数学方法。我们相信，通过刻苦钻研，努力学习，掌握了上述方法后，必将大大提高建立数学模型和运用计算机解决实际问题的综合能力，也一定会在全国建模大赛中取得较好成绩。,培训方式 由于培训内容众多且有相当难度，而培训时间又较短，所以如果学生不事先认真预习相关内容，那么在课堂上不可能完全听懂老师所讲内容，大部分学生会云山雾罩，一头雾水，从而使得培训效果大打折扣。综上，强烈建议学生课前认真、反复研读授课 PPT和其它相关资料，然后带着,疑问和兴趣再听老师讲解，这样才能保证培训效果。培训以讲解为主，适当穿插提问、讨论。学生普遍疑惑的问题，也可当堂咨询老师。,上机软件与上机练习 众所周知，在建模竞赛中，能否熟练使用相关数学软件是能否取得好成绩的关键之一。因此，数学软件的培训应该是建模培训的重要内容。建模中常用的数学软件有Matlab,Spss,Lingo,Maple等。由于培训时间有限，在课堂上只能简单介绍Spss和Lingo，而Matlab 和Maple只,能在上机过程中穿插介绍。几乎所有的数学方法最终都要用数学软件和程序实现，所以上机练习是与课堂讲授同等重要的培训内容。上机训练的主要目的和内容是通过练习，掌握实现各类数学方法的数学软件和程序。上机所用数学软件和程序主要由教师提供，学生只要会用即可。但需要提醒同,学们注意的是，现成的软件和程序不可能解决建模中的所有问题。通过上机训练，掌握一些基本编程和计算技能(如用matlab做数据处理、画图,用maple做简单的解析计算)，对于参加建模竞赛是绝对必要的。上机前，辅导老师会布置上机练习，提供相关软件或程序，讲解关键步骤和程序语句。上机过程中，辅导老师负责解答学生的疑难问题。,第一讲 插值与拟合,插值与拟合属数值分析中函数逼近内容。在数学建模竞赛中，插值与拟合是一种常用的数据分析手段，被公认为建模中的十大算法之一。本节首先通过具体问题引出插值问题与拟合问题，然后简要介绍Matlab中的插值和拟合的相关命令，最后给出两个应用插值和拟合的建模实例。,本节要求学生：(1)理解插值问题和拟合问题；在实际中会正确地判断、选择插值或拟合方法。(2)了解高次插值的Runge 现象及避免方法。(3)熟悉Matlab中一维插值(interp1)、二维插值(interp2)、散乱点插值(griddata)及相关命令(surf,mesh,meshgrid,contour)。(4)熟悉Matlab中多项式拟合(polyfit)、,最小二乘曲线拟合(lsqcurvefit)命令。(5)掌握Matlab编程的基本知识与技能,如数组及运算、调用，循环与控制语句，绘图相关命令，函数(m文件)的定义和调用等。,一、插值问题与拟合问题,引例 矿井中某处的瓦斯浓度 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米已超出所给数据范围，用插值函数外推插值区间外的数据会产生较大的误差。,解决第二个问题的常用方法是，根据地面到井下 500 处的数据求出瓦斯浓度与地面到井下距离x之间的近似函数关系f(x),由f(x)求井下600米处的瓦斯浓度。插值函数过已知点，拟合函数不一定过已知点。通常,插值主要用于求函数值，而拟合的主要目的是求函数关系。当然，某些问题既可以用插值也可以用拟合。,二、高次插值中的Runge现象,通常选用多项式作为插值函数。在研究插值问题的初期，所有人都认为插值多项式的次数越高，插值精度越高。Runge 通过对一个例子的研究发现，上述结论仅仅在插值多项式的次数不超过 七时成立；插值多项式的次数超过七时，插值多项式会出现严重的振荡现象，称之为Runge现象。,例1，节点，求插值多项式。用Maple(Matlab高次插值功能较弱)可方便地求出120次插值多项式，通过图形观察插值效果，见Maple程序演示。,Runge现象,避免 Runge 现象的常用方法是：将插值区间分成若干小区间，在小区间内用低次(二次，三次)插值，即分段低次插值，如样条函数插值。,样条插值结果,三、Matlab插值,Maple 和 Matlab 等数学软件可方便地进行一维和二维多项式插值和样条插值，其中Matlab的二维插值功能较强。Maple中的插值和样条插值命令分别为interp和spline。例如，interp(1,3,4,7,3,5,4,9,x);spline(1,3,4,7,3,5,4,9,x,cubic)。下面介绍Matlab中的插值命令。,1.一维插值 一维插值的典型命令是interp1,其基本格式为yi=interp1(x,y,xi,method)。x,y为插值点，xi,yi为被插值点和插值结果，x,y和xi,yi通常为向量；method表示插值方法：nearest最邻近插值,linear线性插值,spline三次样条插值,cubic立方插值，缺省为线性插值。,例2 在一天24小时内，从零点开始每间隔2小时测得的环境温度数据分别为12，9，9，10，18，24，28，27，25，20，18，15，13，推测中午1点温度，并做出24小时温度曲线图。,Matlab程序 x=0:2:24;y=12 991018 2428272520 18 15 13;x1=13;y1=interp1(x,y,x1,spline)xi=0:1/3600:24;yi=interp1(x,y,xi,spline);plot(x,y,*,xi,yi)运行结果：y1=27.87,请理解、掌握程序中的每个语句，并改变插值方法，观察图形变化。,例3 已知飞机下轮廓线上数据如下，分别画出高次插值(Lagrange)、分段线性插值、样条插值的飞机下轮廓线。,Matlab程序function planex0=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;y1=lagrange(x0,y0,x);y2=interp1(x0,y0,x);y3=interp1(x0,y0,x,spline);subplot(3,1,1),plot(x0,y0,k+,x,y1,r)gridtitle(lagrange)subplot(3,1,2)plot(x0,y0,k+,x,y2,r)gridtitle(piecewise linear)subplot(3,1,3)plot(x0,y0,k+,x,y3,r),gridtitle(spline)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,例3的程序较例2复杂，现说明如下：(1)由于Matlab中没有Lagrange高次插值功能(Maple有)，所以程序中单独编写了高次插值函数lagrange，然后调用；(2)程序中使用了子图形(subplot)、添加网格(grid)和标题(title)、循环和条件语句等。请各位通过上机理解、掌握上述命令,特别是函数的定义及调用。,2.二维插值 二维插值的典型命令是interp2,其基本格式为zi=interp2(x,y,z,xi,yi,method)。二维插值命令的理解和使用较复杂。x,y,z为插值点，z可以理解为被插值函数在(x,y)处的值；xi,yi为被插值点,zi为输出的插值结果，可理解为插值函数在(xi,yi)处的值；x,y为向量，xi,yi为向量或矩阵，而z和zi则为矩阵。,method表示插值方法：nearest最邻近插值,linear双线性插值,spline双三次样条插值,cubic双立方插值,缺省为双线性插值。,例4 测得平板表面3*5网格点处的温度分别为：82 81 80 82 84 79 63 61 65 81 84 84 82 85 86 做出平板表面的温度分布曲面z=f(x,y)的图形及等温线，并求出温度最高和最低点。,Matlab程序x=1:5;y=1:3;temps=82 81 80 82 84;79 63 61 65 81;84 84 82 85 86;figure(1);mesh(x,y,temps)xi=1:0.2:5;yi=1:0.2:3;,zi=interp2(x,y,temps,xi,yi,cubic)figure(2);mesh(xi,yi,zi)figure(3);contour(xi,yi,zi,20,r);i,j=find(zi=min(min(zi);x=xi(j),y=yi(i),zmin=zi(i,j)i,j=find(zi=max(max(zi);x=xi(j),y=yi(i),zmax=zi(i,j),上述程序较复杂，说明如下：(1)interp2中的xi为行向量,而yi为列向量，其实xi和yi行列不同即可。(2)mesh和contour是二维插值及绘图中的常用命令：plot3(空间曲线),mesh(空间曲面)和surf(空间曲面)是3 维作图中的常用命令。mesh 和surf的区别是：mesh画的是曲面网格图，而surf画的是曲面表面图。,contour(x,y,z,n)的功能是作出由若干点(x,y,z)插值而成的曲面的n条等高线。用meshc和surfc也可在曲面下方画出等高线。meshz和surfz是画垂帘图。(3)程序的最后部分为求最高(低)点，请各位通过上机揣摩min,max特别是find的功能。(4)将程序中 figure 语句去除，通过观察结果，体味figure的作用。,例5 在某山区测得一些地点的高程如下表。平面区域为 0 x5600,0y4800试用Matalb中的最邻近插值、双线性插值和双三次插值 3种方法作出该山区的地貌图和等高线图，并求出最高和最低点。,Matlab程序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)xlabel(X),ylabel(Y),zlabel(Z)xi,yi=meshgrid(0:50:5600,0:50:4800);figure(2)z1i=interp2(x,y,z,xi,yi,nearest);surfc(xi,yi,z1i)xlabel(X),ylabel(Y),zlabel(Z)figure(3)z2i=interp2(x,y,z,xi,yi);,surfc(xi,yi,z2i)xlabel(X),ylabel(Y),zlabel(Z)figure(4)z3i=interp2(x,y,z,xi,yi,cubic);surfc(xi,yi,z3i)xlabel(X),ylabel(Y),zlabel(Z)figure(5)subplot(1,3,1),contour(xi,yi,z1i,10,r);subplot(1,3,2),contour(xi,yi,z2i,10,r);subplot(1,3,3),contour(xi,yi,z3i,10,r);,程序中用“xi,yi=meshgrid(0:50:5600,0:50:4800);”生成网格点(xi,yi)，其作用相当于“xi=0:50:5600;yi=0:50:4800;”，可以根据习惯选择其中之一。但 meshgrid(x,y)生成的xi,yi为同维矩阵，xi的行均为x，而yi的列均为y。请各位根据x,y=meshgrid(1:5,2:6)的运行结果理解meshgrid的作用。请补充程序中求最高、高低点部分。,网线图,最邻近点插值,双线性插值,双三次插值,等高线,3.散乱点插值,前面讨论的插值问题的插值点(x,y)均为网格点。当插值点(x,y)为散乱点时，可用命令griddata(x,y,z,xi,yi,method)进行二维插值。例6 在某海域测得一些点(x,y)处的水深z由下表给出，船的吃水深度为5英尺，在矩形区域(75,200)*(-50,150)里的哪些地方船要避免进入。,Matlab程序clearx=129 140 103.5 88 185.5 195 105.5 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;xi,yi=meshgrid(75:0.5:200,-70:0.5:150);zi=griddata(x,y,z,xi,yi,cubic);,figure(1);meshz(xi,yi,zi)xlabel(X),ylabel(Y),zlabel(Z);figure(2),contour(xi,yi,zi,-5-5,b);gridhold onplot(x,y,+)xlabel(X),ylabel(Y),海底地形图,等高线图,四、离散数据的拟合,1.拟合问题及求解 引例 在某化学反应中,根据实验测得生成物浓度y与时间x 的关系如下表，求浓度y与时间x 的对应函数关系y=f(x)，并据此求出反应速度曲线。,显然，从理论上讲 y=f(x)是客观存在的，但在实际中仅由离散数据(xi,yi)是不可能得出y=f(x)的精确表达式的，只能寻找 y=f(x)的一个近似表达式，这种问题称为离散数据的曲线拟合问题。,曲线拟合需解决如下两个问题：(1)线型的选择；(2)中参数的计算。通常主要根据专业知识和散点图确定的线型。参数的计算多采用最小二乘法。Maple、Matlab和Origin等软件都可进行函数拟合，以Origin最为专业。,2.Matlab拟合,Matlab中的拟合命令主要有两个：(1)多项式拟合P=polyfit(x,y,n)x,y为已知数据点，n为指定的多项式次数，P为拟合多项式的所有系数。可用polyval(P,x0)求拟合多项式P在x0处的值。,例7 多项式拟合。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;P=polyfit(x,y,2);z=polyval(P,x);plot(x,y,r*,x,z,b)text(0.5,6,n=2);改变多项式次数，观察拟合结果。,(2)最小二乘拟合P=lsqcurvefit(f,x0,x,y)f为选定的线型，x0为f表达式中参数的初始值，x和y为已知数据点，P为输出项，包括f中的参数、拟合残差等。注意：参数初始值x0对迭代计算是否收敛及收敛速度有影响。另外，不同版本的Matlab中的拟合命令也有差异。,例8 对引例中的数据进行拟合。首先做出数据的散点图。,根据散点图，可选择双指数线型 y=A1*exp(-x/t1)+A2*exp(-x/t2)+y0 下面给出两组计算结果：优初始值结果：z=-2.80 9.31-9.21 1.47 11.16 res=0.0491 劣初始值结果：z=-4.81 2.35-4.81 2.35 10.42 res=0.4166,Matlab程序function myfitx0=1,1,1,1,1;x=1:16;y=4 6.4 8 8.8 9.22 9.5 9.7 9.86 10 10.2 10.32 10.42 10.5 10.55 10.58 10.6;z,res=lsqcurvefit(myfun,x0,x,y)xi=1:0.1:16;yi=myfun(z,xi);plot(x,y,*,xi,yi)function y=myfun(a,x)y=a(1)*exp(-x/a(2)+a(3)*exp(-x/a(4)+a(5);,程序中对拟合函数编制了函数myfun，主函数lsqcurvefit(myfun,x0,x,y)语句对其进行了调用。称为函数句柄，当函数做为参数被调用时，函数名前要加此符号。Matlab拟合过程比较复杂，计算结果受初始值影响较大，精度不高。建议通过上机练习初步掌握Origin 拟合方法。,Origin拟合结果,Equation y=A1*exp(-x/t1)+A2*exp(-x/t2)+y0Adj.R-Square0.99864ValueStandard ErrorBy011.324160.71521BA1-2.766120.23031Bt111.239718.0397BA2-9.351890.51014Bt21.511270.16765,Origin拟合图形,五、上机练习,1.学习、理解、掌握例2例8中的程序，并上机调试。2.在1-12的11小时内，每隔1小时测量一次温度，测得的温度依次为：5,8,9,15,25,29,31,30,22,25,27,24。(1)试估计每隔1/10小时的温度值，并做出温度变化图形。(2)推测t=13.5时的温度。,3.在某山区测得一些地点的高程如下表。平面区域为 1200 x4000,1200y3600试用最邻近插值、双线性插值和双三次插值三种插值方法作出该山区的地貌图和等高线图，并求出最高和最低点。,4.在某海域测得一些点(x,y)处的水深z由下表给出。(1)做出(x,y)的散点图；(2)做出海底地形地貌图；(3)做出伪色彩图(pcolor,shading)；(4)做出彩色等高线图。,(1.486,3.059,0.1);(2.121,4.041,0.1);(2.570,3.959,0.1);(3.439,4.396,0.1);(4.505,3.012,0.1);(3.402,1.604,0.1);(2.570,2.065,0.1);(2.150,1.970,0.1);(1.794,3.059,0.2);(2.121,3.615,0.2);(2.570,3.473,0.2);(3.421,4.160,0.2);(4.271,3.036,0.2);(3.411,1.876,0.2);(2.561,2.562,0.2);(2.179,2.420,0.2);(2.757,3.024,0.3);(3.439,3.970,0.3);(4.084,3.036,0.3);(3.402,2.077,0.3);(2.879,3.036,0.4);(3.421,3.793,0.4);(3.953,3.036,0.4);(3.402,2.219,0.4);(3.000,3.047,0.5);(3.430,3.639,0.5);(3.822,3.012,0.5);(3.411,2.385,0.5);(3.103,3.012,0.6);(3.430,3.462,0.6);(3.710,3.036,0.6);(3.402,2.562,0.6);(3.224,3.047,0.7);(3.411,3.260,0.7);(3.542,3.024,0.7);(3.393,2.763,0.7)。,5.在某化学反应中,测得生成物浓度y与时间x的关系如下表，试分别用多项式插值、非线性最小二乘插值和Origin求浓度 y与时间x 的对应函数关系y=f(x),并据此求出反应速度曲线。,1991年美国大学生数学建模A题,水塔流量的估计,某社区自来水由一个高12.2米,直径17.4米圆柱形水塔提供。当水位降至8.2米时,水泵自动启动加水；水位升高到10.8 米时，水泵停止工作；一般水泵每天工作两次。下表给出了某天不同时间水位数据，其中有三次观察时水泵正在供水，无水位记录。试建立适当的数学模型，计算任意时刻的水流速度，估计一天的用水量和水泵的工作功率。,1.问题分析 题目中给出是水位与时间的关系，而问题的关键是水流速度。由于水塔为圆柱形，水流速度即为水位与时间的函数的导数，所以可采用如下思路获取水流速度：根据给出的水位与时间的对应数据，计算出每个小区间内水流的平均速度，将此速度近似为该小区间中点的瞬时水流速度，从而得到水流速度在各小区间中点的近似值。,对水流速度在各小区间中点的近似值进行拟合或插值，即可得出水流速度与时间的近似关系式。最后对此式积分即得每天的总用水量。至于水泵的功率，可用上述计算结果根据机械原理获得。,2.具体计算方法与结果 可用多项式拟合、样条插值、分段样条插值求取水流速度。拟合时最好根据题给数据选择适当的多项式次数。对于不同方法的计算结果，可根据题给数据进行检验。计算过程及结果见Maple程序演示。本题曾做为校内竞赛试题，冠军张小波给出了极为出色的答卷。,张小波 安徽科艾网络技术有限公司董事长兼总经理，淮南市信息协会副会长，淮南市“云海战略”专家组成员，安徽理工大学计算机学院外聘教师。,多项式拟合结果,用水总量1261，水泵功率1.47。,样条插值结果,用水总量1259，水泵功率1.48。,3、分段样条插值 用水总量1261，水泵功率1.48。,分段样条插值结果,用水总量1261，水泵功率1.48。,2011年全国大学生数学建模A题,城市表层土壤重金属污染分析,随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。,按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、5类区，不同的区域环境受人类活动影响的程度不同。现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距 1公里左右的网格子区域，按照每平方公里 一个采样点对表层土进行取样、编号，并,用GPS记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照 2 公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。附件 1 列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度。,现要求你们通过数学建模来完成以下任务：(1)给出8 种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。(2)通过数据分析，说明重金属污染的主要原因。(3)分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。,(4)分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？,以下是问题(1)的评阅要点：可以用插值的方法获得各重金属污染物的空间分布，再参考由背景值确定的阈值，定量分析城区各区域的污染程度。由于空间数据是不规则的，较好的方法是用散乱数据插值，比如Kriging插值、Shepard插值等。但应明确所使用的方法，并作出分析，不能只简单套用软件。,以下是Matlab插值程序：data=xlsread(data);x=data(:,1);y=data(:,2);z=data(:,4);xi,yi=meshgrid(0:30:28000,0:30:18000);zi=griddata(x,y,z,xi,yi,cubic);,figure(1);mesh(xi,yi,zi)xlabel(X),ylabel(Y),zlabel(Z)figure(2);contour(xi,yi,zi,20,20);grid,Cr空间分布,等高线,Pb空间分布,等高线,