书中Matlab源程序.docx
《书中Matlab源程序.docx》由会员分享，可在线阅读，更多相关《书中Matlab源程序.docx（57页珍藏版）》请在冰豆网上搜索。
书中Matlab源程序
书中Matlab源程序
第1章绪论
【例1-1】有一名学生，期末有5门功课要考试，可用的复习时间有18小时。
假定这五门课程分别是数学、英语、计算机基础、画法几何和专业概论。
如果不复习直接参加考试，这五门功课预期的考试成绩分别为65分、60分、70分、60分和65分。
复习以1小时为一单元，每增加1小时复习时间，各门功课考试成绩就有可能提高，每复习1小时各门功课考试成绩提高的分数分别为3分、4分、5分、4分和6分。
问如何安排各门功课的复习时间可使平均成绩不低于80分，并且数学和英语成绩分别不低于70分和75分。
解：
设分配在数学、英语、计算机基础、画法几何和专业概论这五门功课的复习时间分别为
，则可列出如下的目标函数和限制条件为：
本例具体程序如下：
%li_1_1
f=[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.6667
3.7500
5.0000
0.0000
5.8333
fval=
16.2500
【例1-2】某工厂要生产两种规格的电冰箱，分别用Ⅰ和Ⅱ表示。
生产电冰箱需要两种原材料A和B，另外需设备C。
生产两种电冰箱所需原材料、设备台时、资源供给量及两种产品可获得的利润如表1-1所示。
问工厂应分别生产Ⅰ、Ⅱ型电冰箱多台，才能使工厂获利最多？
表1.1资源需求与限制
资源
Ⅰ
Ⅱ
资源限制
设备
1
1
1200台时
原料A
2
1
1800千克
原料B
0
1
1000千克
单位产品获利
220元
250元
求最大收益
产品Ⅰ用原料限制
800千克
解：
设生产Ⅰ、Ⅱ两种产品的数量分别为
。
则可获得的最大收益为
Matlab求解程序如下：
%li_1_2
clc;
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;
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'）;
holdon
z=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】绘制下面函数的曲线。
解：
应用plot（）函数绘制该函数曲线的程序如下：
%li_1_3
f=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】用图形表示如下优化模型，并求解。
解：
该绘制目标函数曲面、约束函数曲线及求解程序如下：
（1）绘制目标函数曲面的程序
%li_1_4_1
functionli_1_4_1（）
clc;
clearall;
closeall;
n=20;
x1=linspace（0,2,n）;
x2=linspace（0,6,n）;
[xm1,xm2]=meshgrid（x1,x2）;
fori=1:
n
forj=1:
n
xx=[xm1（i,j）,xm2（i,j）];
zm（i,j）=fun_obj（xx）;
end
end
figure
（1）
meshc（xm1,xm2,zm）
xlabel（'x1'）;
ylabel（'x2'）;
zlabel（'zm'）
（2）绘制约束函数曲线及求解的程序
%li_1_4_2
functionli_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:
n
forj=1:
n
xx=[xm1（i,j）,xm2（i,j）];
zm（i,j）=fun_obj（xx）;
end
end
figure
（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）;
holdon
plot（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];
第2章优化设计的数学基础
【例2-11】用图解法分析下面面优化问题的凸性，并求其最小值。
（2-28）
解：
根据所给目标函数和约束函数函数，编制如下程序：
functionkt_test1_plot1
clc;
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）;
holdon
r=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）;
holdon
xx=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]）
xlabel（'x'）;ylabel（'y'）;zlabel（'z'）;
holdoff
figure
（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）;
holdoff
xlabel（'x'）;ylabel（'y'）;
axis（[0605]）;
【例2-12】分析式（2-29）和式（2-30）所示非线性有约束最小值问题。
（2-29）
（2-30）
解：
首先绘制目标函数和约束函数曲面和曲线。
以式（2-29）为例，绘制函数曲面和曲线的程序如下：
functionkt_test1_plot2
clc;
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）;
mesh（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）;
holdon
axis（[-2,2,-2,2,-15,12]）
xlabel（'x'）;ylabel（'y'）;zlabel（'z'）;
holdoff
figure
（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）;
holdon
yy=[-2:
0.1:
2]
plot（zeros（length（yy））,yy,'--'）;
holdon
plot（yy,zeros（length（yy））,'--'）
holdoff
xlabel（'x'）;ylabel（'y'）;
axis（[-22-22]）;
第3章线性规划
【例3-1】用单纯形法求解下面的线性规划问题。
解：
（1）先用MATLAB线性规划函数求解，其计算程序如下：
%ch31_li1
clc;
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_sigma>0
forj=1:
n
sigma（j）=cf（j）-sum（cb（:
）.*A（:
j））;
end
max_sigma=max（sigma）;
if（max_sigma>0）
pvj=find（sigma==max_sigma）;
theta=b./A（:
pvj）;
min_theta=min（theta）;
max_sigma
min_theta
pvi=find（theta==min_theta）;
cb（pvi）=cf（pvj）;
indexb（pvi）=pvj;
pvi
pvj
cb
cf
fori=1:
m
if（i~=pvi）
forj=1:
n
AA（i,j）=A（i,j）-A（i,pvj）*A（pvi,j）/A（pvi,pvj）;
end
bb（i）=b（i）-A（i,pvj）*b（pvi）/A（pvi,pvj）;
else
AA（i,:
）=A（i,:
）/A（pvi,pvj）;
bb（i）=b（i）/A（pvi,pvj）;
end
end
end
A=AA;b=bb';
end
s=1:
n;
x=zeros（n,1）;
fori=1:
m
k=find（s==indexb（i））;
if（k~=0）
x（k）=b（i）;
end
end
fmax=cf*x;
【例3-4】用单纯形法Matlab程序求解例1-2。
解：
为方便起见，重新列出所给问题线性规划的标准形：
编写用户程序。
为目标函数中变量系数行向量；
为初选基变量行向量；indexb1为初始基变量下标索引；
为等式约束方程系数矩阵；
为等式约束方程右端项。
应户程序为：
functionlinear_pro_max_test1
clc
clearall;
cf=[2202500000];
cb=[cf（3）,cf（4）,cf（5）,cf（6） ];
indexb1=[3456];
A=[111000
210100
100010
010001];
b=[1200;1800;800;1000];
[x,fmax]=linear_pro_max（cf,cb,A,b,indexb1）
计算结果：
x=[2001000]。
3.5改进的单纯形法的Matlab程序及实例
根据修正的单纯形法计算步骤，应用Matlab语言编写计算程序，程序清单如下。
function[xfmax]=linear_promax_rev（cf,cb,A,b,indexa1,indexb1）
cf1=cf;
n=length（cf）;
m=length（cb）;
indexa=indexa1;
indexb=indexb1;
x0=b;
forj=1:
n
A1（j）={A（:
j）};
end
e0=[A1{m:
n}];
e0inv=inv（e0）;
xe0=e0inv*x0;
r0=1;
kk=0;
whiler0>0
fprintf（'==========================================================='）
kk=kk+1
indexa,indexb
cb
e0inv
fori=1:
n-m
k=indexa（i）;
cc=A1{k};
r（i）=cf（k）-cb*e0inv*cc;
end
r0=min（r）;
max_sigma=max（r）;
pvj=find（r==max_sigma）;
rr=e0inv*A1{pvj};
theta=x0./rr;
min_theta=min（theta）;
pvi=find（theta==min_theta）;
temp=cb（pvi）;
cb（pvi）=cf（pvj）;
cf（pvj）=temp;
indexa（pvj）=indexb（pvi）;
indexb（pvi）=pvj;
e=A1{pvj};
e1=e0;
fori=1:
m
e1（i,pvi）=e（i）;
end
e1inv=inv（e1）;
xe1=e1inv*b;
e0=e1;
e0inv=e1inv;
x0=xe1;
max_sigma,min_theta,pvi,pvj
e1
e1inv
fprintf（'==========================================================='）
end
s=1:
n;
x=zeros（n,1）;
fori=1:
m
k=find（s==indexb（i））;
if（k~=0）
x（k）=xe1（i）;
end
end
fmax=cf1*x;
【例3-6】应用单纯形法Matlab程序计算例3-5。
解：
将线性规划问题写成标准型，编写用户程序。
为目标函数中变量系数行向量；
为初选基变量行向量；indexa1为初始非基变量下标索引；indexb1为初始基变量下标索引；
为等式约束方程系数矩阵；
等式约束方程右端项。
用户程序如下：
functionlinear_promax_rev_test2
clc;
clearall;
cf=[50100000];
A=[11100
21010
01001];
b=[300;400;250];
cb=[cf（3）,cf（4）cf（5）];
indexa1=[12];
indexb1=[345];
[xfmax]=linear_promax_rev（cf,cb,A,b,indexa1,indexb1）
计算结果为：
x=[502500500],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,th）
%用进退法搜索三个点，使中点函数值最小；输出步长，函数值，设计变量值
%xk0:
初始点
%th:
步长
t1=0;t2=th;
xk1=xk0;xk2=xk1+t2*dir0;
x0=xk1;
f1=feval（f,x0）;
x0=xk2;
f2=feval（f,x0）;
iff2t3=t2+th;
xk3=xk1+t3*dir0;
x0=xk3;
f3=feval（f,x0）;
else
th=-2*th;t3=t1;f3=f1;t1=t2;f1=f2;t2=t3;f2=f3;
t3=th;
xk3=xk1+t3*dir0;
x0=xk3;
f3=feval（f,x0）;
end
ii=0;
whilef2>f3
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;
ifth<0
opt_step=[t3t2t1];
xx=[xx3xx2xx1];
fo=[f3f2f1];
else
opt_step=[t1t2t3];
xx=[xx1xx2xx3];
fo=[f1f2f3];
end
end
【例4-1】用进退法计算函数
的单峰区间，初始点
。
解：
用户程序如下：
%range_search_test1
clc;
f=inline（'2+x^2','x'）;
xk0=2;th=0.5;dir0=1;
[opt_step,fo,xx]=opt_range_serach1（f,xk0,dir0,th）
opt_step=（此句应删除）
结果为：
opt_step=[-3-2-1]
xo=[-101]
fo=[323]
1.黄金分割法的计算框图
图4.5是黄金分割法的计算框图。
图中
为给定的任意小的精度，
为区间缩短的次数。
2.MATLAB程序
%golden_search
function[xo,fo]=golden_search（f,a,b,r,TolX,TolFun,k）
kk=1;
whilekk>0
h=b-a;rh=r*h;c=b-rh;d=a+rh;
fc=feval（f,c）;fd=feval（f,d）;
ifk<=0||abs（h）iffc<=fd
xo=c;fo=fc;
kk=0;
else
xo=d;fo=fd;
kk=0;
end
ifk==0;fprintf（'达到计算次数'）;kk=0;end
else
iffcb=d;k=k-1;
else
a=c;k=k-1;
end
end
end
【例4-2】计算目标函数
在区间
内的极小点。
解：
用户程序如下：
%golden_s_test1.m
functiongolden_s_test1
clc;
clearall;
gs_fun=inline（'2+x^2','x'）;
a=-2;b=2;r=（sqrt（5）-1）/2;TolX=1e-7;TolFun=1e-7;MaxIter=50;
[xo,fo]=fminbnd（gs_fun,a,b）
[xo,fo]=golden_search（gs_fun,a,b,r,TolX,TolFun,MaxIter）
计算结果：
xo=1.4142e-0082;fo=2
2.二次函数插值法的迭代过程与程序框图
MATLAB程序如下：
function[opt_step,fo,xo]=opt_step_quad2（f,xk0,dir0,th,TolX,TolFun,MaxIter）
%opt_step_quad.m
[t012,fo,xx]=opt_range_serach1（f,xk0,dir0,th）;
%searchfortheoptimumstepcorrespondingtominimumf（x）byquadraticapproximationmethod
k=MaxIter;
whilek>0
k=k-1;
t0=t012
（1）;t1=t012
（2）;t2=t012（3）;
x0=xk0+t
升级会员 