常用方法MATLAB求解.ppt

《常用方法MATLAB求解.ppt》由会员分享，可在线阅读，更多相关《常用方法MATLAB求解.ppt（66页珍藏版）》请在三一办公上搜索。

1、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.5 1.0 1.5 2.0 2.5 3.0;,y=1.75 2.45 3.

2、81 4.80 7.00 8.60;,scatter(x,y,5),结果见图,4,由图可看出可用二次多项式拟合。,再输入命令：,p=polyfit(x,y,2),p=0.5614 0.8287 1.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)=

3、yi,i=0,1,n.,解决此问题的拉格朗日插值多项式公式如下,其中Li(x)为n次多项式：,拉格朗日(Lagrange)插值,9,拉格朗日(Lagrange)插值,特别地:,两点一次(线性)插值多项式:,三点二次(抛物)插值多项式:,10,分段线性插值,计算量与n无关;n越大，误差越小.,11,比分段线性插值更光滑。,在数学上，光滑程度的定量描述是：函数(曲线)的k阶导数存在且连续，则称该曲线具有k阶光滑性。光滑性的阶次越高，则越光滑。是否存在较低次的分段多项式达到较高阶光滑性的方法？三次样条插值就是一个很好的例子。,三次样条插值,12,三次样条插值,g(x)为被插值函数。,13,2、一维插

4、值的MATLAB软件命令：,已知离散点上的数据集,求得一解析函数连接自变量相邻的两个点，并求得两点,间的数值，这一过程叫插值。,MATLAB在一维插值函数interp1中，提供了四种插值方法选择：线性插值、三次样条插值、立方插值和最近邻点插值。interp1的本格式为：,yi=interp1(x，y，xi，method),其中x,y分别表示数据点的横、纵坐标向量，x 必须单调，xi为需要插值的横坐标数据（或数组），xi不能超出x的范围，而method为可选参数，有四种选择：,nearest：最邻近插值,linear：线性插值；,14,spline：三次样条插值；,cubic：立方插值。,缺省时

5、：分段线性插值。,例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=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=10.2000 30.0000 30.9000 24.9000,15,T=interp1(hours,temps,3.2 6.5 7.1 11.7,spline)%三次样条

6、插值,T=9.6734 30.0427 31.1755 25.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分别给出了针对插值基点为网格节点的插值函数及针对插值基点为散乱节点的插值函数调用格式

7、。,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 三次

8、插值,缺省时,是线性插值,19,例3：测得平板表面35网格点处的温度分别为：,82 81 80 82 84 79 63 61 65 81 84 84 82 85 86,试求在平板表面坐标为（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=82 81 80 82 84;79 63 61 65 81;84 84 82 85 86;,20,mesh(x0,y0,temps),结果见图5.6,分别用线性性插

9、值和三次样条插值求已知点的温度。,输入命令：,21,t=interp2(x0,y0,temps,1.5 2 2.5 3.5,1.5 1.6 2 4.5,liner),t=76.2500 70.2000 62.0000 NaN,T=interp2(x0,y0,temps,1.5 2 2.5 3.5,1.5 1.6 2 4.5,spline),T=71.4531 65.5200 60.9688 188.8906,（2）以平滑数据,在x、y方向上每隔0.2个单位的地方进行插值画出线性和三次样条插值的温度分布曲面图.,输入以下命令得温度的线性插值曲面图:,x=1:0.2:5;,y=1:0.2:3;,z

10、=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,1200=y=3600)试

11、作出该山区的地貌图和等高线图，并对几种插值方法进行比较。,25,此例将对最近邻点插值、线性插值方法和双三次插值方法的插值效果进行比较。,输入命令：,x0=0:400:5600;,y0=0:400:4800;,z0=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

12、 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 1

13、550 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

14、 950 820 690 540 380 300 210;.1350 1370 1390 1400 1410 960 940 880 800 690 570 430 290 210 150;,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,每隔

15、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是待求

16、函数值所给定网格点的的横、纵坐标向量。,而method为可选参数，有四种选择：,nearest 最邻近插值,linear 线性插值,cubic 三次插值,30,v4 Matlab提供的插值方法,缺省时,线性插值,例5 在某海域测得一些点(x,y)处的水深z由下表5.3给出，船的吃水深度为5英尺，在矩形区域（75，200）*（-50，150）里的哪些地方船要避免进入。,31,(1)作出测量点的分布图,MATLAB命令：,x0=129 140 103.5 88 185.5 195 105.5 157.5 107.5 77 81 162 162 117.5;,y0=7.5 141.5 23 147

17、22.5 137.5 85.5-6.5-81 3 56.5-66.5 84-33.5;,plot(x0,y0,+);,结果见图5.13,32,（2）作出三次插值海底地形图,输入命令：,x0=129 140 103.5 88 185.5 195 105.5 157.5 107.5 77 81 162 162 117.5;,y0=7.5 141.5 23 147 22.5 137.5 85.5-6.5-81 3 56.5-66.5 84-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,cub

18、ic);,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软件求解一些优化问

19、题，包括求解线性规划和非线性规划等问题。,线性规划是运筹学的一个重要分支，它起源于工业生产组织管理的决策问题。在数学上它用来确定多变量线性函数在变量满足线性约束条件下的最优值；随着计算机的发展，出现了如单纯形法等有效算法，它在工农业、军事、交通运输、决策管理与规划等领域中有广泛的应用。,37,问题一：,任务分配问题：某车间有甲、乙两台机床，可用于加工三种工件。假定这两台车床的可工作时间分别为800和900，三种工件的数量分别为400、600和500，且已知车床甲加工单位数量三种工件所需的时间和加工费分别为0.4、1.1、1和13、9、10，车床乙加工单位数量三种工件所需的时间和加工费分别为0.

20、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元。为使总检验费用最省，该工厂应聘一级、二级检验员各几名？,解

21、 设需要一级和二级检验员的人数分别为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,模型求解

22、,最优解：x14=x21=x32=x43=1,其它变量为0;成绩为253.2(秒)=413”2,MIN 66.8x11+75.6x12+87x13+58.6x14+67.4x51+71 x52+83.8x53+62.4x54SUBJECT TO x11+x12+x13+x14=1 x41+x42+x43+x44=1 x11+x21+x31+x41+x51=1 x14+x24+x34+x44+x54=1END INT 20,输入LINDO求解,甲 自由泳、乙 蝶泳、丙 仰泳、丁 蛙泳.,44,丁蛙泳c43=69.675.2，戊自由泳c54=62.4 57.5,方案是否调整？,乙 蝶泳、丙 仰泳、

23、丁 蛙泳、戊 自由泳,最优解：x21=x32=x43=x51=1,成绩为417”7,c43,c54 的新数据重新输入模型，用LINDO求解,指派(Assignment)问题：每项任务有且只有一人承担，每人只能承担一项，效益不同，怎样分派使总效益最大.,讨论,45,为了选修课程门数最少，应学习哪些课程？,问题四 选课策略,要求至少选两门数学课、三门运筹学课和两门计算机课,选修课程最少，且学分尽量多，应学习哪些课程？,46,0-1规划模型,决策变量,目标函数,xi=1 选修课号i 的课程（xi=0 不选）,选修课程总数最少,约束条件,最少2门数学课，,3门运筹学课，,2门计算机课。,47,先修课程

24、要求,最优解：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,

25、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=cX s.t AXb,命令：x=linprog（c，A，b）,2、模型：min z=cX,命令：x=li

26、nprog(c,A,b,Aeq,beq),如果没有不等式：,存在，则令A=,b=.,55,3、模型：min z=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.01 0.01 0.01 0.03 0.03 0.03;0.02 0 0 0.05 0 0;0 0.02 0 0 0.05 0;0 0 0.03 0

27、0 0.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.5000 0.5000 3.0000 0.0000 0.0000 0.0000,fval=-2.5000e+004,58,例2,s.t,0,i=1,2 j=1,2,3.,解 用命令3，编写M文件xxgh2.m如下：,c=10 5 6 4 8 12;,A=;,b=;,Aeq=1 1 1 0 0 0;0 0 0 1 1 1;1 0 0 1 0 0;0 1 0 0 1

28、 0;0 0 1 0 0 1;,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.0000 20.0000 40.0000 50.0000 50.0000 0.0000,fval=940.0000,60,例3 用MATLAB求问题一：,编写M文件xxgh3.m如下:,c=13 9 10 11 12 8;,A=0.4 1.1 1 0 0 0；0 0 0 0.5 1.2 1.3;,b=800;900;,Aeq=1 0 0 1 0 0；0 1 0 0 1 0；0 0

29、1 0 0 1;,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.0000 600.0000 0.0000 400.0000 0.0000 500.0000,fval=1.3800e+004,即在甲机床上加工600个工件2,在乙机床上加工400个工件1、500个工件3，可在满足条件的情况下使总加工费最小为13800。,62,4、一般非线性规划求解,线性规划标准型为：,min F(X)s.t AX=b,G(X),Ceq(X)=0 VLB,X,VUB,其中X为n维变元向

30、量，G(X)与Ceq(X)均为非线性函数组成的向量，其它变量的含义与线性规划、二次规划中相同.用MATLAB求解上述问题，基本步骤分三步：,1.首先建立M文件fun.m,定义目标函数F（X）:,function f=fun(X);,f=F(X);,63,2、若约束条件中有非线性约束:G(X),则建立M文件nonlcon.m定义函数G(X)与Ceq(X):,或Ceq(X)=0,function G,Ceq=nonlcon(X),G=.,Ceq=.,3.建立主程序.非线性规划求解的函数是fmincon,命令的基本格式如下：,x,fval=fmincon(fun,X0,A,b,Aeq,beq,VLB

31、,VUB,nonlcon),64,x,fval=fmincon(fun,X0,A,b,Aeq,beq,VLB,VUB,nonlcon),说明：fan为目标函数的表达式字符串；如果没有不等式：AXb存在，则令A=,b=；如果没有等式约束：AeqX=beq存在，则令Aeq=,beq=；如果没有不等式：vlbxvub约束，则令vlb=,ulb=；nonlcon的作用是通过接收的向量x来计算非线性不等约束G(x)0和等式约束Ceq(x)=0分别在x处的估计G和Ceq,通过指定函数柄来使用；函数可能会出现局部最优解，这与初值x0的选取有关。,65,例,1）先建立M文件 fun4.m,定义目标函数:,fu

32、nction f=fun4(x);,f=exp(x(1)*(4*x(1)2+2*x(2)2+4*x(1)*x(2)+2*x(2)+1);,2）再建立M文件mycon.m定义非线性约束：,function g,ceq=mycon(x),g=1.5+x(1)*x(2)-x(1)-x(2);-x(1)*x(2)-10;,ceq=;,66,3）主程序youh3.m为:,x0=-1;1;,A=;b=;,Aeq=1 1;beq=0;,vlb=;vub=;,x,fval=fmincon(fun4,x0,A,b,Aeq,beq,vlb,vub,mycon),3.运算结果为：,x=-1.2250 1.2250,fval=1.8951,