常用方法MATLAB求解

1,数学建模竞赛,2,一、曲线拟合及MATLAB软件求解,已知离散点上的数据集,求得一解析函数y=f(x)使y=f(x)在原离散点,接近给定,曲线拟合是最小二乘法曲线拟合，拟合结果可使误差的,上尽可能,的值，这一过程叫曲线拟合。最常用的,平方和最小，即找出使,最小的f(x).,几种常见的数学方法及软件求解,3,格式：p=polyfit(x,y,n).,说明：求出已知数据x,y的n次拟合多项式f(x)的系数p,x必须是单调的。,例1已知某函数的离散值如表,求二次拟合多项式.,先画函数离散点的图形,输入命令：,x=0.51.01.52.02.53.0;,y=1.752.453.814.807.008.60;,scatter(x,y,5),结果见图,4,由图可看出可用二次多项式拟合。,再输入命令：,p=polyfit(x,y,2),p=0.56140.82871.1560,即二次拟合多项式为,5,画出离散点及拟合曲线：,输入命令：,x1=0.5:0.05:3.0;,y1=polyval(p,x1);,plot(x,y,*r,x1,y1,-b),结果见图5.4,6,1、一维插值的定义,二、一维插值,7,8,称为拉格朗日插值基函数。,已知函数f(x)在n+1个点x0,x1,xn处的函数值为y0,y1,yn。求一n次多项式函数Pn(x)，使其满足：Pn(xi)=yi,i=0,1,n.,解决此问题的拉格朗日插值多项式公式如下,其中Li(x)为n次多项式：,拉格朗日(Lagrange)插值,9,拉格朗日(Lagrange)插值,特别地:,两点一次(线性)插值多项式:,三点二次(抛物)插值多项式:,10,分段线性插值,计算量与n无关;n越大，误差越小.,11,比分段线性插值更光滑。,在数学上，光滑程度的定量描述是：函数(曲线)的k阶导数存在且连续，则称该曲线具有k阶光滑性。光滑性的阶次越高，则越光滑。是否存在较低次的分段多项式达到较高阶光滑性的方法？三次样条插值就是一个很好的例子。,三次样条插值,12,三次样条插值,g(x)为被插值函数。,13,2、一维插值的MATLAB软件命令：,已知离散点上的数据集,求得一解析函数连接自变量相邻的两个点，并求得两点,间的数值，这一过程叫插值。,MATLAB在一维插值函数interp1中，提供了四种插值方法选择：线性插值、三次样条插值、立方插值和最近邻点插值。interp1的本格式为：,yi=interp1(x，y，xi，method),其中x,y分别表示数据点的横、纵坐标向量，x必须单调，xi为需要插值的横坐标数据（或数组），xi不能超出x的范围，而method为可选参数，有四种选择：,nearest：最邻近插值,linear：线性插值；,14,spline：三次样条插值；,cubic：立方插值。,缺省时：分段线性插值。,例2在1-12的11小时内，每隔1小时测量一次温度，测得的温度依次为：5，8，9，15，25，29，31，30，22，25，27，24。试估计在3.2,6.5,7.1,11.7小时的温度值。,解输入命令：,hours=1:12;,temps=589152529313022252724;,t=interp1(hours,temps,3.26.57.111.7)%线性插值,t=10.200030.000030.900024.9000,15,T=interp1(hours,temps,3.26.57.111.7,spline)%三次样条插值,T=9.673430.042731.175525.3820,比较发现，两种结果有差异，这是因为插值是一个估计或猜测的过程。,两种插值的画图如下；,输入命令：,t0=1:0.1:12;,T0=interp1(hours,temps,t0,spline);,plot(hours,temps,+,t0,T0,hours,temps,r:),xlabel(时间);,ylabel(温度),16,gtext(线性插值),gtext(三次样条插值),结果见图5.5,17,三、二维插值,对二维插值问题，MATLAB分别给出了针对插值基点为网格节点的插值函数及针对插值基点为散乱节点的插值函数调用格式。,1、用MATLAB作网格节点数据的插值,对上述问题，MATLAB提供了二维插值函数interp2,其基本格式为：,18,z=interp2(x0,y0,z0,x,y,method),其中x0,y0是自变量。X0,y0的分量值必须是单调递增的。X0和y0分别是m维和n维向量，分别表示已知数据点的横、纵坐标向量，z0是m*n维矩阵，标明相应于所给数据网格点的函数值。向量x,y是待求函数值所给定网格点的的横、纵坐标向量，x,y的值分别不能超出x0,y0的范围。,而method为可选参数，有四种选择：,nearest最邻近插值,linear线性插值,spline三次样条插值,cubic三次插值,缺省时,是线性插值,19,例3：测得平板表面35网格点处的温度分别为：,828180828479636165818484828586,试求在平板表面坐标为（1.5,1.5）,（2,1.6）,（2.5,2）（3.5,4.5）处的温度，并作平板表面的温度分布曲面z=f(x,y)的图形，,（1）先在三维坐标画出原始数据，画出粗糙的温度分布曲图.,输入以下命令：,x0=1:5;,y0=1:3;,temps=8281808284;7963616581;8484828586;,20,mesh(x0,y0,temps),结果见图5.6,分别用线性性插值和三次样条插值求已知点的温度。,输入命令：,21,t=interp2(x0,y0,temps,1.522.53.5,1.51.624.5,liner),t=76.250070.200062.0000NaN,T=interp2(x0,y0,temps,1.522.53.5,1.51.624.5,spline),T=71.453165.520060.9688188.8906,（2）以平滑数据,在x、y方向上每隔0.2个单位的地方进行插值画出线性和三次样条插值的温度分布曲面图.,输入以下命令得温度的线性插值曲面图:,x=1:0.2:5;,y=1:0.2:3;,z=interp2(x0,y0,temps,x,y,linear);,mesh(x,y,z),22,xlabel(x轴);,ylabel(y轴);,zlabel(z轴),title(线性插值曲面图),结果见图5.7,23,再输入以下命令得温度的三次样条插值曲面图:,z=interp2(x0,y0,temps,x,y,spline);,mesh(x,y,z),xlabel(x轴);,ylabel(y轴);,zlabel(z轴),title(三次样条插值曲面图),结果见图5.7,24,例4山区地貌：,在某山区测得一些地点的高程如下表。平面区域为1200=x=4000,1200x0=0:400:5600;,y0=0:400:4800;,z0=370470550600670690670620580450400300100150250;.510620730800850870850780720650500200300350320;.650760880970102010501020830900700300500550480350;.740880108011301250128012301040900500700780750650550;.830980118013201450142014001300700900850840380780750;.88010601230139015001500140090011001060950870900930950;.9101090127015001200110013501450120011501010880100010501100;.9501190137015001200110015501600155013801070900105011501200;.143014301460150015501600155016001600160015501500150015501550;.1420143014501480150015501510143013001200980850750550500;.138014101430145014701320128012001080940780620460370350;.13701390141014301440114011101050950820690540380300210;.13501370139014001410960940880800690570430290210150;,26,meshz(x0,y0,z0),xlabel(x轴);,ylabel(y轴);,zlabel(z轴),title(原始数据山区表面图),结果见图5.8,27,（1）线性插值,输入命令：,z2=interp2(x0,y0,z0,x,y,linear);,surfc(x,y,z2),xlabel(x轴);,ylabel(y轴);,zlabel(z轴),title(线性插值表面图),结果见图5.10,每隔50加密网格，分别作线性插值、三次样条插值，作出插值后的表面图：,28,（2）三次样条插值,输入命令：,z3=interp2(x0,y0,z0,x,y,spline);,surfc(x,y,z3),xlabel(x轴);,ylabel(y轴);,zlabel(z轴),title(三次样条插值表面图),结果见图5.11,29,2、用MATLAB作散乱节点数据的插值,对上述问题，MATLAB提供了二维插值函数griddata,其基本格式为：,z=griddata（x0，y0，z0，x，y，method）,其中X0、y0、z0均是n维向量，分别表示已知数据点的横、纵坐标和竖坐标向量。向量x,y是待求函数值所给定网格点的的横、纵坐标向量。,而method为可选参数，有四种选择：,nearest最邻近插值,linear线性插值,cubic三次插值,30,v4Matlab提供的插值方法,缺省时,线性插值,例5在某海域测得一些点(x,y)处的水深z由下表5.3给出，船的吃水深度为5英尺，在矩形区域（75，200）*（-50，150）里的哪些地方船要避免进入。,31,(1)作出测量点的分布图,MATLAB命令：,x0=129140103.588185.5195105.5157.5107.57781162162117.5;,y0=7.5141.52314722.5137.585.5-6.5-81356.5-66.584-33.5;,plot(x0,y0,+);,结果见图5.13,32,（2）作出三次插值海底地形图,输入命令：,x0=129140103.588185.5195105.5157.5107.57781162162117.5;,y0=7.5141.52314722.5137.585.5-6.5-81356.5-66.584-33.5;,z0=-4-8-6-8-6-8-8-9-9-8-8-9-4-9;,y=-70:0.5:150;,z=griddata(x0,y0,z0,x,y,cubic);,meshz(x,y,z);,xlabel(x轴);,ylabel(y轴);,zlabel(z轴),title(三次插值海底地形图),x=75:0.5:200;,33,结果见图,34,（3）作出水深低于5英尺的海域范围,输入命令：,contour(x,y,z,-5,-5,-k);,grid,xlabel(x轴);,ylabel(y轴);,title(船不能进入区域),结果见图,35,从图可以看出，船应避免进入危险区域,（-10，110）。,（110，190）,36,四、MATLAB软件在优化中的应用,在研究与解决具体问题中，经常遇到有关优化问题，下面介绍用MATLAB软件求解一些优化问题，包括求解线性规划和非线性规划等问题。,线性规划是运筹学的一个重要分支，它起源于工业生产组织管理的决策问题。在数学上它用来确定多变量线性函数在变量满足线性约束条件下的最优值；随着计算机的发展，出现了如单纯形法等有效算法，它在工农业、军事、交通运输、决策管理与规划等领域中有广泛的应用。,37,问题一：,任务分配问题：某车间有甲、乙两台机床，可用于加工三种工件。假定这两台车床的可工作时间分别为800和900，三种工件的数量分别为400、600和500，且已知车床甲加工单位数量三种工件所需的时间和加工费分别为0.4、1.1、1和13、9、10，车床乙加工单位数量三种工件所需的时间和加工费分别为0.5、1.2、1.3和11、12、8。问怎样分配车床的加工任务，才能既满足加工工件的要求，又使加工费用最低？,38,解设在甲车床上加工工件1、2、3的数量分别为x1、x2、x3，在乙车床上加工工件1、2、3的数量分别为x4、x5、x6。,可建立以下线性规划模型：,39,问题二：某厂每日8小时的产量不低于1800件。为了进行质量控制，计划聘请两种不同水平的检验员。一级检验员的标准为：速度25件/小时，正确率98%，计时工资4元/小时；二级检验员的标准为：速度15件/小时，正确率95%，计时工资3元/小时。检验员每错检一次，工厂要损失2元。为使总检验费用最省，该工厂应聘一级、二级检验员各几名？,解设需要一级和二级检验员的人数分别为x1、x2人,则应付检验员的工资为：,因检验员错检而造成的损失为：,故目标函数为：,40,约束条件为：,线性规划模型：,41,丁的蛙泳成绩退步到115”2；戊的自由泳成绩进步到57”5,组成接力队的方案是否应该调整？,如何选拔队员组成4100米混合泳接力队?,问题三混合泳接力队的选拔,5名候选人的百米成绩,穷举法：组成接力队的方案共有5!=120种。,42,目标函数,若选择队员i参加泳姿j的比赛，记xij=1,否则记xij=0,0-1规划模型,cij(秒)队员i第j种泳姿的百米成绩,约束条件,每人最多入选泳姿之一,每种泳姿有且只有1人,43,模型求解,最优解：x14=x21=x32=x43=1,其它变量为0;成绩为253.2(秒)=413”2,MIN66.8x11+75.6x12+87x13+58.6x14+67.4x51+71x52+83.8x53+62.4x54SUBJECTTOx11+x12+x13+x14=1x41+x42+x43+x44=1x11+x21+x31+x41+x51=1x14+x24+x34+x44+x54=1ENDINT20,输入LINDO求解,甲自由泳、乙蝶泳、丙仰泳、丁蛙泳.,44,丁蛙泳c43=69.675.2，戊自由泳c54=62.457.5,方案是否调整？,乙蝶泳、丙仰泳、丁蛙泳、戊自由泳,最优解：x21=x32=x43=x51=1,成绩为417”7,c43,c54的新数据重新输入模型，用LINDO求解,指派(Assignment)问题：每项任务有且只有一人承担，每人只能承担一项，效益不同，怎样分派使总效益最大.,讨论,45,为了选修课程门数最少，应学习哪些课程？,问题四选课策略,要求至少选两门数学课、三门运筹学课和两门计算机课,选修课程最少，且学分尽量多，应学习哪些课程？,46,0-1规划模型,决策变量,目标函数,xi=1选修课号i的课程（xi=0不选）,选修课程总数最少,约束条件,最少2门数学课，,3门运筹学课，,2门计算机课。,47,先修课程要求,最优解：x1=x2=x3=x6=x7=x9=1,其它为0；6门课程，总学分21,0-1规划模型,约束条件,x3=1必有x1=x2=1,模型求解（LINDO）,48,学分最多,多目标优化的处理方法：化成单目标优化。,两目标(多目标)规划,讨论：选修课程最少，学分尽量多，应学习哪些课程？,课程最少,以学分最多为目标，不管课程多少。,以课程最少为目标，不管学分多少。,49,多目标规划,在课程最少的前提下以学分最多为目标。,最优解：x1=x2=x3=x5=x7=x9=1,其它为0；总学分由21增至22。,注意：最优解不唯一！,LINDO无法告诉优化问题的解是否唯一。,可将x9=1易为x6=1,50,多目标规划,对学分数和课程数加权形成一个目标，例如三七开。,最优解：x1=x2=x3=x4=x5=x6=x7=x9=1，其它为0；总学分28。,51,讨论与思考,最优解与1=0，2=1的结果相同学分最多,多目标规划,最优解与1=1，2=0的结果相同课程最少,52,一般线性规划问题的数学表达式：,max(min)f=,s.t,.,用矩阵向量符号表示：,max(min)f=cX,53,其中,max(min)f=cX,54,用MATLAB优化工具箱解线性规划的命令如下：,1、模型：minz=cXs.tAXb,命令：x=linprog（c，A，b）,2、模型：minz=cX,命令：x=linprog(c,A,b,Aeq,beq),如果没有不等式：,存在，则令A=,b=.,55,3、模型：minz=cX,命令：x=linprog(c,A,b,Aeq,beq,vlb,vub),如果没有等式约束：,则令Aeq=,beq=.,存在，,4、命令：x,fval=linprog(),返回最优解及处的目标函数值fval.,56,解用命令3，编写M文件xxgh1.m如下：,c=-0.4-0.28-0.32-0.72-0.64-0.6;,A=0.010.010.010.030.030.03;0.02000.0500;00.02000.050;000.03000.08;,b=850;700;100;900;,Aeq=;beq=;,vlb=0;0;0;0;0;0;vub=;,57,x,fval=linprog(c,A,b,Aeq,beq,vlb,vub),结果：,x=1.0e+004*3.50000.50003.00000.00000.00000.0000,fval=-2.5000e+004,58,例2,s.t,0,i=1,2j=1,2,3.,解用命令3，编写M文件xxgh2.m如下：,c=10564812;,A=;,b=;,Aeq=111000;000111;100100;010010;001001;,59,beq=60;100;50;70;40;,vlb=0;0;0;0;0;0;vub=;,x,fval=linprog(c,A,b,Aeq,beq,vlb,vub),结果：,x=0.000020.000040.000050.000050.00000.0000,fval=940.0000,60,例3用MATLAB求问题一：,编写M文件xxgh3.m如下:,c=1391011128;,A=0.41.11000；0000.51.21.3;,b=800;900;,Aeq=100100；010010；001001;,61,beq=400;600;500;,vlb=0;0;0;0;0;0;,vub=;,x,fval=linprog(c,A,b,Aeq,beq,vlb,vub),结果:,x=0.0000600.00000.0000400.00000.0000500.0000,fval=1.3800e+004,即在甲机床上加工600个工件2,在乙机床上加工400个工件1、500个工件3，可在满足条件的情况下使总加工费最小为13800。,62,4、一般非线性规划求解,线性规划标准型为：,minF(X)s.tAX=b,G(X),Ceq(X)=0VLB,X,VUB,其中X为n维变元