书中Matlab源程序.docx

1、书中Matlab源程序书中Matlab源程序第1 章 绪论【例1-1】有一名学生，期末有5门功课要考试，可用的复习时间有18小时。假定这五门课程分别是数学、英语、计算机基础、画法几何和专业概论。如果不复习直接参加考试，这五门功课预期的考试成绩分别为65分、60分、70分、60分和65分。复习以1小时为一单元，每增加1小时复习时间，各门功课考试成绩就有可能提高，每复习1小时各门功课考试成绩提高的分数分别为3分、4分、5分、4分和6分。问如何安排各门功课的复习时间可使平均成绩不低于80分，并且数学和英语成绩分别不低于70分和75分。解：设分配在数学、英语、计算机基础、画法几何和专业概论这五门功课的

2、复习时间分别为，则可列出如下的目标函数和限制条件为： 本例具体程序如下：%li_1_1f=1 1 1 1 1;A=1 1 1 1 1; -3 -4 -5 -4 -6; -3 0 0 0 0; 0 -4 0 0 0; 3 0 0 0 0; 0 4 0 0 0; 0 0 5 0 0; 0 0 0 4 0; 0 0 0 0 6;b=18;-80;-5;-15;35;40;30;40;35;lb=zeros(6,1)x,fval=linprog(f,A,b,lb)计算结果为：x = 1.6667 3.7500 5.0000 0.0000 5.8333fval = 16.2500 【例1-2】某工厂要生

3、产两种规格的电冰箱，分别用和表示。生产电冰箱需要两种原材料A和B，另外需设备C。生产两种电冰箱所需原材料、设备台时、资源供给量及两种产品可获得的利润如表1-1所示。问工厂应分别生产、型电冰箱多台，才能使工厂获利最多？表1.1 资源需求与限制资源资源限制设备111200台时原料A211800千克原料B011000千克单位产品获利220元250元求最大收益产品用原料限制800千克解：设生产、两种产品的数量分别为。则可获得的最大收益为Matlab求解程序如下：%li_1_2clc;close all;f=-220 250;A=1 1;2 1;1 0;0 1;b=1200;1800;800;1000;

4、xl=0 0;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;plot(x1,x21,x1,x22,0:1:1000,1000,800,0:1:1500,x1,x23,r)axis(0,1400,0,2000)xlabel(x1);ylabel(x2);hold onz=200*xm1+250*xm2;C,h=contour(xm1,xm2,z);text_handle = clabel(C,h);se

5、t(text_handle,BackgroundColor,1 1 .6,Edgecolor,.7 .7 .7);hold off【例1-3】绘制下面函数的曲线。解：应用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.1leq Theta leq 2pi)ylabel(2sin(Theta)+

6、ln(Theta);title(Plot of 2sin(Theta)+ln(Theta)text(pi/4,sin(-pi/4),leftarrow 2sin(Theta)+ln(Theta),HorizontalAlignment,left)legend(-)grid on【例1-4】用图形表示如下优化模型，并求解。 解：该绘制目标函数曲面、约束函数曲线及求解程序如下：（1）绘制目标函数曲面的程序%li_1_4_1function li_1_4_1()clc;clear all;close all;n=20;x1=linspace(0,2,n);x2=linspace(0,6,n);xm1

7、,xm2=meshgrid(x1,x2);for i=1:nfor j=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_2function li_1_4_2()clc;clear all;close all;x0=1,1;x,fval,exitflag,output=fmincon(fun_obj,x0, , , , , , ,fun_cons)n=20;x1=linspace

8、(0,6,n);x2=linspace(0,10,n);xm1,xm2=meshgrid(x1,x2);for i=1:nfor j=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)plot(x1,y1,x1,y2);hold onplot(x1(1:8),y3(1:8),-r)hold

9、onc,h=contour(xm1,xm2,zm,20);clabel(c,h);xlabel(x1);ylabel(x2);function f=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;第2章 优化设计的数学基础【例2-11】 用图解法分析下面面优化问题的凸性，并求其最小值。 （2-28）解：根据所给目标函数和约束函数函数，编制如下程序：function kt_test1_plot1clc;cle

10、ar all;close all;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);hold onr=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);hold onxx=linspace(-1,6,30);yy=linspace(-2,5,30);xxm,zzm=meshgrid(0*xx,0,6);yym=mesh

11、grid(yy,-10,0);mesh(xxm,yym,zzm);hold onyym,zzm=meshgrid(0*yy,0,6);xxm=meshgrid(xx,-10,0);mesh(xxm,yym,zzm);axis(-1,6,-3,5,-2,12)xlabel(x);ylabel(y);zlabel(z);hold offfigure(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);hold onc,h=contour(xm,ym,z

12、m3,10);text=clabel(c);hold offxlabel(x);ylabel(y);axis(0 6 0 5 );【例2-12】分析式（2-29）和式（2-30）所示非线性有约束最小值问题。 （2-29） （2-30）解：首先绘制目标函数和约束函数曲面和曲线。以式（2-29）为例，绘制函数曲面和曲线的程序如下：function kt_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);figu

13、re(1);mesh(xm,ym,zm1);hold onxm,ym=meshgrid(linspace(-2,2,30),linspace(-2,2,30);zm2=(-(xm-1).2-(ym-2).2+10);mesh(xm,ym,zm2);hold onaxis(-2,2,-2,2,-15,12)xlabel(x);ylabel(y);zlabel(z);hold offfigure(2)x=linspace(-1,1,30);y1=sqrt(1-x.2);y2=-sqrt(1-x.2);plot(x,y1,k,x,y2,-r);hold onc,h=contour(xm,ym,zm2

14、,20);text=clabel(c);hold onyy=-2:0.1:2plot(zeros(length(yy),yy,-);hold onplot(yy,zeros(length(yy),-)hold offxlabel(x);ylabel(y);axis(-2 2 -2 2 );第 3 章 线性规划【例3-1】用单纯形法求解下面的线性规划问题。解：（1）先用MATLAB线性规划函数求解，其计算程序如下：%ch31_li1clc;close all;f=-2 1 -2;A=1 1 2;1 3 -1;b=5 3;xl=0 0 0;x,fval=linprog(f,A,b, , ,xl)3

15、.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);while max_sigma0for j=1:n sigma(j)=cf(j)-sum(cb(:).*A(:,j);endmax_sigma=max(sigma);if(max_sigma0)pvj=find(sigma=max_sigma);theta=b./A(:,pvj);min_theta

16、=min(theta);max_sigmamin_thetapvi=find(theta=min_theta);cb(pvi)=cf(pvj);indexb(pvi)=pvj;pvipvjcbcffor i=1:m if(i=pvi)for j=1:nAA(i,j)=A(i,j)-A(i,pvj)*A(pvi,j)/A(pvi,pvj);endbb(i)=b(i)-A(i,pvj)*b(pvi)/A(pvi,pvj);elseAA(i,:)=A(i,:)/A(pvi,pvj);bb(i)=b(i)/A(pvi,pvj);endendendA=AA;b=bb;ends=1:n;x=zeros(n

17、,1);for i=1:mk=find(s=indexb(i);if(k=0)x(k)=b(i);endendfmax=cf*x;【例3-4】用单纯形法Matlab程序求解例1-2。解：为方便起见，重新列出所给问题线性规划的标准形：编写用户程序。为目标函数中变量系数行向量；为初选基变量行向量；indexb1为初始基变量下标索引；为等式约束方程系数矩阵；为等式约束方程右端项。应户程序为：function linear_pro_max_test1clcclear all;cf=220 250 0 0 0 0;cb=cf(3),cf(4),cf(5),cf(6);indexb1=3 4 5 6;A=

18、1 1 1 0 0 0 2 1 0 1 0 0 1 0 0 0 1 00 1 0 0 0 1;b=1200;1800;800;1000;x,fmax=linear_pro_max(cf,cb,A,b,indexb1)计算结果：x=200 1000。3.5改进的单纯形法的Matlab程序及实例 根据修正的单纯形法计算步骤，应用Matlab语言编写计算程序，程序清单如下。function x fmax=linear_promax_rev(cf,cb,A,b,indexa1,indexb1)cf1=cf;n=length(cf);m=length(cb);indexa=indexa1;indexb=

19、indexb1;x0=b;for j=1:nA1(j)=A(:,j);ende0=A1m:n;e0inv=inv(e0);xe0=e0inv*x0;r0=1;kk=0;while r00 fprintf(=)kk=kk+1indexa,indexbcbe0invfor i=1:n-mk=indexa(i);cc=A1k;r(i)=cf(k)-cb*e0inv*cc;endr0=min(r);max_sigma=max(r);pvj=find(r=max_sigma);rr=e0inv*A1pvj;theta=x0./rr;min_theta=min(theta);pvi=find(theta=

20、min_theta);temp=cb(pvi);cb(pvi)=cf(pvj);cf(pvj)=temp;indexa(pvj)=indexb(pvi);indexb(pvi)=pvj;e=A1pvj;e1=e0;for i=1:me1(i,pvi)=e(i);ende1inv=inv(e1);xe1=e1inv*b;e0=e1;e0inv=e1inv;x0=xe1;max_sigma,min_theta,pvi,pvje1e1invfprintf(=)ends=1:n;x=zeros(n,1);for i=1:mk=find(s=indexb(i);if(k=0)x(k)=xe1(i);en

21、dendfmax=cf1*x;【例3-6】 应用单纯形法Matlab程序计算例3-5。解：将线性规划问题写成标准型，编写用户程序。为目标函数中变量系数行向量；为初选基变量行向量；indexa1为初始非基变量下标索引；indexb1为初始基变量下标索引；为等式约束方程系数矩阵；等式约束方程右端项。用户程序如下：function linear_promax_rev_test2clc;clear all;cf=50 100 0 0 0;A=1 1 1 0 0 2 1 0 1 0 0 1 0 0 1;b=300;400;250;cb=cf(3),cf(4) cf(5);indexa1=1 2;inde

22、xb1=3 4 5;x fmax=linear_promax_rev(cf,cb,A,b,indexa1,indexb1)计算结果为：x=50 250 0 50 0,fmax=27500。第 4 章 一维搜索方法2. 进退法的MATLAB程序该程序是按照寻找最佳步长编写的，其数学模型为式（4-2）。程序说明如下：函数opt_range_serach(xk0,dir0,th)中输入参数；xk0：初始点；dir0：给定的搜索方向；th：步长增量。%opt_range_serach1.m function opt_step,fo,xx = opt_range_serach1(f,xk0,dir0,t

23、h)%用进退法搜索三个点，使中点函数值最小；输出步长，函数值，设计变量值 %xk0:初始点%th:步长t1=0;t2=th; xk1=xk0;xk2=xk1+t2*dir0; x0=xk1; f1=feval(f,x0); x0=xk2; f2=feval(f,x0); if f2f3 t1=t2;f1=f2;t2=t3;f2=f3; t3=t2+th; xk3=xk1+t3*dir0; x0=xk3; f3=feval(f,x0); end xx1=xk1+t1*dir0; xx2=xk1+t2*dir0; xx3=xk1+t3*dir0; if th0h = b - a; rh = r*h; c = b - rh; d = a + rh;fc = feval(f,c); fd = feval(f,d);if k = 0 | abs(h) TolX & abs(fc - fd) TolFun if fc = fdxo = c;fo = fc;kk=0;elsexo = d;fo = fd;kk=0;endif k = 0;fprintf(达到计算次数);kk=0; endelseif fc 0k=k-1;t0=t012(1);t1=t012(2);t2=t012(3);x0 = xk0+t012