书中Matlab源程序.docx
- 文档编号:9009709
- 上传时间:2023-05-16
- 格式:DOCX
- 页数:66
- 大小:173.93KB
书中Matlab源程序.docx
《书中Matlab源程序.docx》由会员分享,可在线阅读,更多相关《书中Matlab源程序.docx(66页珍藏版)》请在冰点文库上搜索。
书中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);
iff2 t3=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 iffc b=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+t012
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- Matlab 源程序