Matlab源程序

1书中Matlab源程序第1章绪论【例1-1】有一名学生，期末有5门功课要考试，可用的复习时间有18小时。假定这五门课程分别是数学、英语、计算机基础、画法几何和专业概论。如果不复习直接参加考试，这五门功课预期的考试成绩分别为65分、60分、70分、60分和65分。复习以1小时为一单元，每增加1小时复习时间，各门功课考试成绩就有可能提高，每复习1小时各门功课考试成绩提高的分数分别为3分、4分、5分、4分和6分。问如何安排各门功课的复习时间可使平均成绩不低于80分，并且数学和英语成绩分别不低于70分和75分。解：设分配在数学、英语、计算机基础、画法几何和专业概论这五门功课的复习时间分别为12345,,,,xxxxx，则可列出如下的目标函数和限制条件为：12345()fxxxxxx1234518xxxxx12345(34546320)/580xxxxx136570x本例具体程序如下：%li_1_1f=[11111];A=[11111;-3-4-5-4-6;-30000;0-4000;30000;04000;00500;00040;00006];b=[18;-80;-5;-15;35;40;30;40;35];lb=zeros(6,1)[x,fval]=linprog(f,A,b,[],[],lb)计算结果为：x=1.66673.75005.00000.00005.8333fval=16.25002【例1-2】某工厂要生产两种规格的电冰箱，分别用Ⅰ和Ⅱ表示。生产电冰箱需要两种原材料A和B，另外需设备C。生产两种电冰箱所需原材料、设备台时、资源供给量及两种产品可获得的利润如表1-1所示。问工厂应分别生产Ⅰ、Ⅱ型电冰箱多台，才能使工厂获利最多？表1.1资源需求与限制资源ⅠⅡ资源限制设备111200台时原料A211800千克原料B011000千克单位产品获利220元250元求最大收益产品Ⅰ用原料限制800千克解：设生产Ⅰ、Ⅱ两种产品的数量分别为12,xx。则可获得的最大收益为212max()220250,fxxxxR12121212..1200218008001000,0stxxxxxxxxMatlab求解程序如下：%li_1_2clc;closeall;f=-[220250];A=[11;21;10;01];b=[1200;1800;800;1000];xl=[00];[x,fval]=linprog(f,A,b,[],[],xl)x1=[0:1800];x2=[0:2000];[xm1,xm2]=meshgrid(x1,x2);x21=1200-x1;x22=1800-2*x1;x23=(-fval-220*x1)/250;3plot(x1,x21,x1,x22,[0:1:1000],1000,800,[0:1:1500],x1,x23,'r')axis([0,1400,0,2000])xlabel('x1');ylabel('x2');holdonz=200*xm1+250*xm2;[C,h]=contour(xm1,xm2,z);text_handle=clabel(C,h);set(text_handle,'BackgroundColor',[11.6],'Edgecolor',[.7.7.7]);holdoff【例1-3】绘制下面函数的曲线。)ln()sin(2)(xxxy解：应用plot()函数绘制该函数曲线的程序如下：%li_1_3f=inline('2*sin(x)+log(x)','x')x=linspace(0.1,2*pi,15);y=feval(f,x);plot(x,y,'-rs','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','g','MarkerSize',10)xlabel('0.1\leq\Theta\leq2\pi')ylabel('2sin(\Theta)+ln(\Theta)');title('Plotof2sin(\Theta)+ln(\Theta)')text(pi/4,sin(-pi/4),'\leftarrow2sin(\Theta)+ln(\Theta)','HorizontalAlignment','left')legend('-')gridon【例1-4】用图形表示如下优化模型，并求解。124)(min2231xxxf016)5()5(25..2221221xxxxtsx解：该绘制目标函数曲面、约束函数曲线及求解程序如下：（1）绘制目标函数曲面的程序%li_1_4_1functionli_1_4_1()clc;clearall;closeall;n=20;4x1=linspace(0,2,n);x2=linspace(0,6,n);[xm1,xm2]=meshgrid(x1,x2);fori=1:nforj=1:nxx=[xm1(i,j),xm2(i,j)];zm(i,j)=fun_obj(xx);endendfigure(1)meshc(xm1,xm2,zm)xlabel('x1');ylabel('x2');zlabel('zm')（2）绘制约束函数曲线及求解的程序%li_1_4_2functionli_1_4_2()clc;clearall;closeall;x0=[1,1];[x,fval,exitflag,output]=fmincon(@fun_obj,x0,[],[],[],[],[],[],@fun_cons)n=20;x1=linspace(0,6,n);x2=linspace(0,10,n);[xm1,xm2]=meshgrid(x1,x2);fori=1:nforj=1:nxx=[xm1(i,j),xm2(i,j)];zm(i,j)=fun_obj(xx);endendfigure(1)f1=inline('sqrt(25-x.^2)','x');f2=inline('sqrt(16-(x-5).^2)+5','x');y1=feval(f1,x1);y2=feval(f2,x1);y3=sqrt((4*x1.^3)-12-fval+0.01)5plot(x1,y1,x1,y2);holdonplot(x1(1:8),y3(1:8),'--r')holdon[c,h]=contour(xm1,xm2,zm,20);clabel(c,h);xlabel('x1');ylabel('x2');functionf=fun_obj(x)f=4*x(1)^3-x(2)^2-12;function[c,ceq]=fun_cons(x)c=[x(1)^2+x(2)^2-10*x(1)-10*x(2)+34-x(1)-x(2)];ceq=[x(1)^2+x(2)^2-25];6第2章优化设计的数学基础【例2-11】用图解法分析下面面优化问题的凸性，并求其最小值。221221122132min()(3)..()40()0()0fxxstgxxgxgxxxxx（2-28）解：根据所给目标函数和约束函数函数，编制如下程序：functionkt_test1_plot1clc;clearall;closeall;x=linspace(0,2.5,30);[xm,zm1]=meshgrid(x,[0,6]);y=4-x.^2;ym=meshgrid(y,[0,20]);mesh(xm,ym,zm1);holdonr=2;t=linspace(0,2*pi,45);x=r*cos(t)+3;y=r*sin(t);[xm,ym]=meshgrid(x,y);zm2=(xm-3).^2+ym.^2;mesh(xm,ym,zm2);holdonxx=linspace(-1,6,30);yy=linspace(-2,5,30);[xxm,zzm]=meshgrid(0*xx,[0,6]);[yym]=meshgrid(yy,[-10,0]);mesh(xxm,yym,zzm);holdon[yym,zzm]=meshgrid(0*yy,[0,6]);[xxm]=meshgrid(xx,[-10,0]);mesh(xxm,yym,zzm);axis([-1,6,-3,5,-2,12])7xlabel('x');ylabel('y');zlabel('z');holdofffigure(2)x=linspace(0,6,30);y=linspace(0,5,30);[xm,ym]=meshgrid(x,y);zm3=(xm-3).^2+ym.^2;y1=4-x.^2;plot(x,y1,'k');holdon[c,h]=contour(xm,ym,zm3,10);text=clabel(c);holdoffxlabel('x');ylabel('y');axis([0605]);【例2-12】分析式（2-29）和式（2-30）所示非线性有约束最小值问题。22122211212min()(1)(10)10..()1000fxxstgxxxxxx（2-29）22122211212min()(1)(10)10..()1000fxxstgxxxxxx（2-30）解：首先绘制目标函数和约束函数曲面和曲线。以式（2-29）为例，绘制函数曲面和曲线的程序如下：functionkt_test1_plot2clc;thita=linspace(0,2*pi,50);r=1;x=r*cos(thita);y=r*sin(thita);[xm,zm1]=meshgrid(x,[-5,10]);[ym]=meshgrid(y,[0,20]);figure(1);8mesh(xm,ym,zm1);holdon[xm,ym]=meshgrid(linspace(-2,2,30),linspace(-2,2,30));zm2=(-(xm-1).^2-(ym-2).^2+10);mesh(xm,ym,zm2);holdonaxis([-2,2,-2,2,-15,12])xlabel('x');ylabel('y');zlabel('z');holdofffigure(2)x=linspace(-1,1,30);y1=sqrt(1-x.^2);y2=-sqrt(1-x.^2);plot(x,y1,'k',x,y2,'-r');holdon[c,h]=contour(xm,ym,zm2,20);text=clabel(c);holdonyy=[-2:0.1:2]plot(zeros(length(yy)),yy,'--');holdonplot(yy,zeros(length(yy)),'--')holdoffxlabel('x');ylabel('y');axis([-22-22]);9第3章线性规划【例3-1】用单纯形法求解下面的线性规划问题。0,,3352..22)(max321321321321xxxxxxxxxtsxxxfx解：（1）先用MATLAB线性规划函数求解，其计算程序如下：%ch31_li1clc;closeall;f=-[21-2];A=[112;13-1];b=[53];xl=[000];[x,fval]=linprog(f,A,b,[],[],xl)3.3单纯形法的Matlab程序及实例程序清单如下：function[x,fmax]=linear_pro_max(cf,cb,A,b,indexb1)n=length(cf);max_sigma=1;m=length(cb);indexb=indexb1;theta=zeros(size(m,1));whilemax_sigma0forj=1:nsigma(j)=cf(j)-sum(cb(:).*A(:,j)