常用方法MATLAB求解

常用方法MATLAB求解1第1页，共66页，2023年，2月20日，星期一一、曲线拟合及MATLAB软件求解

已知离散点上的数据集求得一解析函数y=f(x)使y=f(x)在原离散点接近给定曲线拟合是最小二乘法曲线拟合，拟合结果可使误差的上尽可能的值，这一过程叫曲线拟合。最常用的平方和最小，即找出使最小的f(x).几种常见的数学方法及软件求解2第2页，共66页，2023年，2月20日，星期一格式：p=polyfit(x,y,n).说明：求出已知数据x,y的n次拟合多项式f(x)的系数p,x必须是单调的。例1已知某函数的离散值如表xi0.51.01.52.02.53.0yi1.752.453.814.807.008.65求二次拟合多项式.先画函数离散点的图形输入命令：>>x=[0.51.01.52.02.53.0];>>y=[1.752.453.814.807.008.60];>>scatter(x,y,5)结果见图3第3页，共66页，2023年，2月20日，星期一由图可看出可用二次多项式拟合。再输入命令：>>p=polyfit(x,y,2)p=0.56140.82871.1560即二次拟合多项式为4第4页，共66页，2023年，2月20日，星期一画出离散点及拟合曲线：输入命令：>>x1=0.5:0.05:3.0;>>y1=polyval(p,x1);>>plot(x,y,'*r',x1,y1,'-b')

结果见图5.45第5页，共66页，2023年，2月20日，星期一1、一维插值的定义已知n+1个节点其中互不相同，不妨设求任一插值点处的插值二、一维插值6第6页，共66页，2023年，2月20日，星期一构造一个(相对简单的)函数通过全部节点,即再用计算插值，即7第7页，共66页，2023年，2月20日，星期一称为拉格朗日插值基函数。

已知函数f(x)在n+1个点x0,x1,…,xn处的函数值为y0,y1,…,yn。求一n次多项式函数Pn(x)，使其满足：

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

三次样条插值g(x)为被插值函数。12第12页，共66页，2023年，2月20日，星期一2、一维插值的MATLAB软件命令：已知离散点上的数据集求得一解析函数连接自变量相邻的两个点，并求得两点间的数值，这一过程叫插值。MATLAB在一维插值函数interp1中，提供了四种插值方法选择：线性插值、三次样条插值、立方插值和最近邻点插值。interp1的本格式为：yi=interp1(x，y，xi，'method')其中x,y分别表示数据点的横、纵坐标向量，x必须单调，xi为需要插值的横坐标数据（或数组），xi不能超出x的范围，而method为可选参数，有四种选择：‘nearest’：最邻近插值‘linear’：线性插值；13第13页，共66页，2023年，2月20日，星期一‘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.900014第14页，共66页，2023年，2月20日，星期一>>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('温度')15第15页，共66页，2023年，2月20日，星期一>>gtext('线性插值')>>gtext('三次样条插值')结果见图5.516第16页，共66页，2023年，2月20日，星期一三、二维插值对二维插值问题，MATLAB分别给出了针对插值基点为网格节点的插值函数及针对插值基点为散乱节点的插值函数调用格式。1、

用MATLAB作网格节点数据的插值

对上述问题，MATLAB提供了二维插值函数interp2,其基本格式为：17第17页，共66页，2023年，2月20日，星期一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’三次插值缺省时,是线性插值18第18页，共66页，2023年，2月20日，星期一例3：测得平板表面3×5网格点处的温度分别为：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];19第19页，共66页，2023年，2月20日，星期一mesh(x0,y0,temps)结果见图5.6分别用线性性插值和三次样条插值求已知点的温度。输入命令：20第20页，共66页，2023年，2月20日，星期一>>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)21第21页，共66页，2023年，2月20日，星期一>>xlabel('x轴');>>ylabel('y轴');>>zlabel('z轴')>>title('线性插值曲面图')结果见图5.722第22页，共66页，2023年，2月20日，星期一再输入以下命令得温度的三次样条插值曲面图:>>z=interp2(x0,y0,temps,x',y,'spline');

>>mesh(x,y,z)>>xlabel('x轴');>>ylabel('y轴');>>zlabel('z轴')>>title('三次样条插值曲面图')

结果见图5.723第23页，共66页，2023年，2月20日，星期一例4山区地貌：在某山区测得一些地点的高程如下表。平面区域为1200<=x<=4000,1200<=y<=3600)试作出该山区的地貌图和等高线图，并对几种插值方法进行比较。Yx04008001200160020002400280032003600400044004800520056000400800120016002000240028003200360040004400480037047055060067069067062058045040030010015025051062073080085087085078072065050020030035032065076088097010201050102083090070030050055048035074088010801130125012801230104090050070078075065055083098011801320145014201400130070090085084038078075088010601230139015001500140090011001060950870900930950910109012701500120011001350145012001150101088010001050110095011901370150012001100155016001550138010709001050115012001430143014601500155016001550160016001600155015001500155015501420143014501480150015501510143013001200980850750550500138014101430145014701320128012001080940780620460370350137013901410143014401140111010509508206905403803002101350137013901400141096094088080069057043029021015024第24页，共66页，2023年，2月20日，星期一此例将对最近邻点插值、线性插值方法和双三次插值方法的插值效果进行比较。输入命令：>>x0=0:400:5600;>>y0=0:400:4800;>>z0=[370470550600670690670620580450400300100150250;...510620730800850870850780720650500200300350320;...650760880970102010501020830900700300500550480350;...740880108011301250128012301040900500700780750650550;...830980118013201450142014001300700900850840380780750;...88010601230139015001500140090011001060950870900930950;...9101090127015001200110013501450120011501010880100010501100;...9501190137015001200110015501600155013801070900105011501200;...143014301460150015501600155016001600160015501500150015501550;...1420143014501480150015501510143013001200980850750550500;...138014101430145014701320128012001080940780620460370350;...13701390141014301440114011101050950820690540380300210;...13501370139014001410960940880800690570430290210150];25第25页，共66页，2023年，2月20日，星期一>>meshz(x0,y0,z0)>>xlabel('x轴');>>ylabel('y轴');>>zlabel('z轴')>>title('原始数据山区表面图')结果见图5.826第26页，共66页，2023年，2月20日，星期一（1）线性插值输入命令：>>z2=interp2(x0,y0,z0,x,y','linear');>>surfc(x,y,z2)>>xlabel('x轴');>>ylabel('y轴');>>zlabel('z轴')>>title('线性插值表面图')结果见图5.10每隔50加密网格，分别作线性插值、三次样条插值，作出插值后的表面图：27第27页，共66页，2023年，2月20日，星期一（2）三次样条插值输入命令：>>z3=interp2(x0,y0,z0,x,y','spline');>>surfc(x,y,z3)>>xlabel('x轴');>>ylabel('y轴');>>zlabel('z轴')>>title('三次样条插值表面图')结果见图5.1128第28页，共66页，2023年，2月20日，星期一2、

用MATLAB作散乱节点数据的插值

对上述问题，MATLAB提供了二维插值函数griddata,其基本格式为：z=griddata（x0，y0，z0，x，y，‘method’）其中X0、y0、z0均是n维向量，分别表示已知数据点的横、纵坐标和竖坐标向量。向量x,y是待求函数值所给定网格点的的横、纵坐标向量。而method为可选参数，有四种选择：‘nearest’最邻近插值‘linear’线性插值‘cubic’三次插值29第29页，共66页，2023年，2月20日，星期一'v4'Matlab提供的插值方法缺省时,线性插值例5在某海域测得一些点(x,y)处的水深z由下表5.3给出，船的吃水深度为5英尺，在矩形区域（75，200）*（-50，150）里的哪些地方船要避免进入。X129.0140.0103.588.0185.5195.0105.5157.5107.577.081.0162.0162.0117.5Y7.5141.523.0147.022.5137.585.5－6.5－81.03.056.5－66.584.0－33.5Z4868688998894930第30页，共66页，2023年，2月20日，星期一(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.1331第31页，共66页，2023年，2月20日，星期一（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;32第32页，共66页，2023年，2月20日，星期一结果见图33第33页，共66页，2023年，2月20日，星期一（3）作出水深低于5英尺的海域范围输入命令：>>contour(x,y,z,[-5,-5],’-k’);>>grid,>>xlabel('x轴');>>ylabel('y轴');>>title('船不能进入区域')结果见图34第34页，共66页，2023年，2月20日，星期一从图可以看出，船应避免进入危险区域（-10，110）。（110，190）35第35页，共66页，2023年，2月20日，星期一四、MATLAB软件在优化中的应用在研究与解决具体问题中，经常遇到有关优化问题，下面介绍用MATLAB软件求解一些优化问题，包括求解线性规划和非线性规划等问题。线性规划是运筹学的一个重要分支，它起源于工业生产组织管理的决策问题。在数学上它用来确定多变量线性函数在变量满足线性约束条件下的最优值；随着计算机的发展，出现了如单纯形法等有效算法，它在工农业、军事、交通运输、决策管理与规划等领域中有广泛的应用。36第36页，共66页，2023年，2月20日，星期一

问题一：

任务分配问题：某车间有甲、乙两台机床，可用于加工三种工件。假定这两台车床的可工作时间分别为800和900，三种工件的数量分别为400、600和500，且已知车床甲加工单位数量三种工件所需的时间和加工费分别为0.4、1.1、1和13、9、10，车床乙加工单位数量三种工件所需的时间和加工费分别为0.5、1.2、1.3和11、12、8。问怎样分配车床的加工任务，才能既满足加工工件的要求，又使加工费用最低？37第37页，共66页，2023年，2月20日，星期一解设在甲车床上加工工件1、2、3的数量分别为x1、x2、x3，在乙车床上加工工件1、2、3的数量分别为x4、x5、x6。可建立以下线性规划模型：38第38页，共66页，2023年，2月20日，星期一问题二：某厂每日8小时的产量不低于1800件。为了进行质量控制，计划聘请两种不同水平的检验员。一级检验员的标准为：速度25件/小时，正确率98%，计时工资4元/小时；二级检验员的标准为：速度15件/小时，正确率95%，计时工资3元/小时。检验员每错检一次，工厂要损失2元。为使总检验费用最省，该工厂应聘一级、二级检验员各几名？解

设需要一级和二级检验员的人数分别为x1、x2人,则应付检验员的工资为：因检验员错检而造成的损失为：故目标函数为：39第39页，共66页，2023年，2月20日，星期一约束条件为：线性规划模型：40第40页，共66页，2023年，2月20日，星期一丁的蛙泳成绩退步到1’15”2；戊的自由泳成绩进步到57”5,组成接力队的方案是否应该调整？如何选拔队员组成4100米混合泳接力队?问题三混合泳接力队的选拔

甲乙丙丁戊蝶泳1’06”857”21’18”1’10”1’07”4仰泳1’15”61’06”1’07”81’14”21’11”蛙泳1’27”1’06”41’24”61’09”61’23”8自由泳58”653”59”457”21’02”45名候选人的百米成绩穷举法：组成接力队的方案共有5!=120种。41第41页，共66页，2023年，2月20日，星期一目标函数若选择队员i参加泳姿j的比赛，记xij=1,否则记xij=0

0-1规划模型

cij(秒)~队员i第j种泳姿的百米成绩约束条件每人最多入选泳姿之一

ciji=1i=2i=3i=4i=5j=166.857.2787067.4j=275.66667.874.271j=38766.484.669.683.8j=458.65359.457.262.4每种泳姿有且只有1人42第42页，共66页，2023年，2月20日，星期一模型求解

最优解：x14=x21=x32=x43=1,其它变量为0;成绩为253.2(秒)=4’13”2MIN66.8x11+75.6x12+87x13+58.6x14+……+67.4x51+71x52+83.8x53+62.4x54SUBJECTTOx11+x12+x13+x14<=1

……x41+x42+x43+x44<=1x11+x21+x31+x41+x51=1

……x14+x24+x34+x44+x54=1ENDINT20

输入LINDO求解

甲乙丙丁戊蝶泳1’06”857”21’18”1’10”1’07”4仰泳1’15”61’06”1’07”81’14”21’11”蛙泳1’27”1’06”41’24”61’09”61’23”8自由泳58”653”59”457”21’02”4甲~自由泳、乙~蝶泳、丙~仰泳、丁~蛙泳.43第43页，共66页，2023年，2月20日，星期一丁蛙泳c43

=69.675.2，戊自由泳c54=62.4

57.5,方案是否调整？乙~蝶泳、丙~仰泳、丁~蛙泳、戊~自由泳最优解：x21=x32=x43=x51=1,成绩为4’17”7c43,c54的新数据重新输入模型，用LINDO求解指派(Assignment)问题：每项任务有且只有一人承担，每人只能承担一项，效益不同，怎样分派使总效益最大.讨论甲~自由泳、乙~蝶泳、丙~仰泳、丁~蛙泳.原方案44第44页，共66页，2023年，2月20日，星期一为了选修课程门数最少，应学习哪些课程？

问题四选课策略要求至少选两门数学课、三门运筹学课和两门计算机课课号课名学分所属类别先修课要求1微积分5数学

2线性代数4数学

3最优化方法4数学；运筹学微积分；线性代数4数据结构3数学；计算机计算机编程5应用统计4数学；运筹学微积分；线性代数6计算机模拟3计算机；运筹学计算机编程7计算机编程2计算机

8预测理论2运筹学应用统计9数学实验3运筹学；计算机微积分；线性代数选修课程最少，且学分尽量多，应学习哪些课程？

45第45页，共66页，2023年，2月20日，星期一0-1规划模型

决策变量

目标函数

xi=1~选修课号i的课程（xi=0~不选）

选修课程总数最少约束条件课号课名所属类别1微积分数学2线性代数数学3最优化方法数学；运筹学4数据结构数学；计算机5应用统计数学；运筹学6计算机模拟计算机；运筹学7计算机编程计算机8预测理论运筹学9数学实验运筹学；计算机最少2门数学课，3门运筹学课，2门计算机课。46第46页，共66页，2023年，2月20日，星期一先修课程要求最优解：

x1=x2=x3=x6=x7=x9=1,其它为0；6门课程，总学分210-1规划模型

约束条件x3=1必有x1=x2=1模型求解（LINDO）课号课名先修课要求1微积分

2线性代数

3最优化方法微积分；线性代数4数据结构计算机编程5应用统计微积分；线性代数6计算机模拟计算机编程7计算机编程

8预测理论应用统计9数学实验微积分；线性代数47第47页，共66页，2023年，2月20日，星期一学分最多多目标优化的处理方法：化成单目标优化。两目标(多目标)规划

讨论：选修课程最少，学分尽量多，应学习哪些课程？课程最少以学分最多为目标，不管课程多少。以课程最少为目标，不管学分多少。最优解如上，6门课程，总学分21。最优解显然是选修所有9门课程。48第48页，共66页，2023年，2月20日，星期一多目标规划

在课程最少的前提下以学分最多为目标。最优解：

x1=x2=x3=x5=x7=x9=1,其它为0；总学分由21增至22。注意：最优解不唯一！课号课名学分1微积分52线性代数43最优化方法44数据结构35应用统计46计算机模拟37计算机编程28预测理论29数学实验3LINDO无法告诉优化问题的解是否唯一。可将x9=1易为x6=1增加约束，以学分最多为目标求解。49第49页，共66页，2023年，2月20日，星期一多目标规划

对学分数和课程数加权形成一个目标，例如三七开。最优解：

x1=x2=x3=x4=x5=x6=x7=x9=1，其它为0；总学分28。课号课名学分1微积分52线性代数43最优化方法44数据结构35应用统计46计算机模拟37计算机编程28预测理论29数学实验350第50页，共66页，2023年，2月20日，星期一讨论与思考最优解与1=0，2=1的结果相同——学分最多多目标规划

最优解与1=1，2=0的结果相同——课程最少51第51页，共66页，2023年，2月20日，星期一一般线性规划问题的数学表达式：

max(min)f=s.t

,

,

……

.用矩阵向量符号表示：max(min)f=cX52第52页，共66页，2023年，2月20日，星期一其中

max(min)f=cX53第53页，共66页，2023年，2月20日，星期一用MATLAB优化工具箱解线性规划的命令如下：1、模型：minz=cXs.tAX≤b命令：x=linprog（c，A，b）

2、模型：minz=cX

命令：x=linprog(c,A,b,Aeq,beq)如果没有不等式：存在，则令A=[],b=[].54第54页，共66页，2023年，2月20日，星期一3、模型：minz=cX

命令：x=linprog(c,A,b,Aeq,beq,vlb,vub)如果没有等式约束：则令Aeq=[],beq=[].存在，4、命令：[x,fval]=linprog(…)返回最优解ｘ及ｘ处的目标函数值fval.55第55页，共66页，2023年，2月20日，星期一解用命令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=[];56第56页，共66页，2023年，2月20日，星期一[x,fval]=linprog(c,A,b,Aeq,beq,vlb,vub)结果：x=1.0e+004*3.50000.50003.00000.00000.00000.0000fval=-2.5000e+00457第57页，共66页，2023年，2月20日，星期一

例2

s.t

0,i=1,2j=1,2,3.解用命令3，编写M文件xxgh2.m如下：c=[10564812];A=[];b=[];Aeq=[111000;000111;100100;010010;001001];58第58页，共66页，2023年，2月20日，星期一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.0000fval=940.000059第59页，共66页，2023年，2月20日，星期一例3用MATLAB求问题一：编写M文件xxgh3.m如下:c=[1391011128];A=[0.41.1