工程一班23号甘超第三次作业.docx
《工程一班23号甘超第三次作业.docx》由会员分享，可在线阅读，更多相关《工程一班23号甘超第三次作业.docx（10页珍藏版）》请在冰豆网上搜索。
工程一班23号甘超第三次作业
第三次作业
P1266-2
一
（1）混合函数法
所用到的牛顿法M文件
function[x,minf]=minNT（f,x0,var,eps）
formatlong;
ifnargin==3
eps=1.0e-6;
end
tol=1;
x0=transpose（x0）;
gradf=jacobian（f,var）;
jacf=jacobian（gradf,var）;
whiletol>eps
v=subs（gradf,var,x0）;
tol=norm（v）;
pv=subs（jacf,var,x0）;
p=-inv（pv）*transpose（v）;
p=double（p）;
x1=x0+p;
x0=x1;
end
x=x1;
minf=subs（f,var,x）;
formatshort;
所用到的混合函数M文件
function[x,minf]=minMixFun6-2（f,g,x0,r0,c,var,eps）
gx0=subs（g,var,x0）;
ifgx0>=0;
else
disp（'³õʼµã±ØÐëÂú×ã²»µÈʽԼÊø!
'）;
x=NaN;
minf=NaN;
return;
end
ifr0<=0
disp（'³õʼÕÏ° Òò×Ó±ØÐë´óÓÚ0!
'）;
x=NaN;
minf=NaN;
return;
end
ifc>=1||c<0
disp（'ËõСϵÊý±ØÐë´óÓÚ0ÇÒСÓÚ1!
'）;
x=NaN;
minf=NaN;
return;
end
ifnargin==7
eps=1.0e-6;
end
FE=0;
fori=1:
length（g）
FE=FE+1/g（i）;
end
FH=transpose（h）*h;
x1=transpose（x0）;
x2=inf;
while1
FF=r0*FE+FH/sqrt（r0）;
SumF=f+FF;
[x2,minf]=minNT（SumF,transpose（x1）,var）;
ifnorm（x2-x1）<=eps
x=x2;
break;
else
r0=c*r0;
x1=x2;
end
end
minf=subs（f,var,x）;
输入
>>symsst;
>>f=s^2+2*t^2;
>>g=[s+t-1;s;t];
>>h=[s+t-1;s;t];
>>[x,minf]=minMixFun（f,g,h,[1,1],2,0.5,[s,t],1.0e-6）
结果
x=
0.3333
0.3333
minf=
0.3333
二
（2）内点罚函数法
所插入的牛顿法M文件
function[x,minf]=minNT（f,x0,var,eps）
formatlong;
ifnargin==3
eps=1.0e-6;
end
tol=1;
x0=transpose（x0）;
gradf=jacobian（f,var）;
jacf=jacobian（gradf,var）;
whiletol>eps
v=subs（gradf,var,x0）;
tol=norm（v）;
pv=subs（jacf,var,x0）;
p=-inv（pv）*transpose（v）;
p=double（p）;
x1=x0+p;
x0=x1;
end
x=x1;
minf=subs（f,var,x）;
formatshort;
内点惩罚函数法M文件
function[x,minf]=minNF（f,x0,g,u,v,var,eps）
formatlong;
ifnargin==6
eps=1.0e-6;
end
k=0;
FE=0;
fori=1:
length（g）
FE=FE+1/g（i）;
end
x1=transpose（x0）;
x2=inf;
while1
FF=u*FE;
SumF=f+FF;
[x2,minf]=minNT（SumF,transpose（x1）,var）;
Bx=subs（FE,var,x2）;
ifu*Bxifnorm（x2-x1）<=eps
x=x2;
break;
else
u=u*v;
x1=x2;
end
else
ifnorm（x2-x1）<=eps
x=x2;
break;
else
u=u*v;
x1=x2;
end
end
end
minf=subs（f,var,x）;
formatshort;
输入
>>clear
>>symst;a=10;b=5;
>>f=a*t;
>>g=[t-b;t];
>>[x,minf]=minNF（f,[8],g,20,0.5,[t],1.0e-6）
输出结果
x=
5.0000
minf=
50.0000
单纯形法求解线性规划问题
单纯形法M文件
function[x,minf]=simpleMthd（A,c,b,baseVector）
sz=size（A）;
nVia=sz
（2）;
n=sz
（1）;
xx=1:
nVia;
nobase=zeros（1,1）;
m=1;
fori=1:
nVia;
if（isempty（find（baseVector==xx（i）,1）））
nobase（m）=i;
m=m+1;
else
;
end
end
bCon=1;
M=0;
whilebCon
nB=A（:
nobase）;
ncb=c（nobase）;
B=A（:
baseVector）;
cb=c（baseVector）;
xb=inv（B）*b;
f=cb*xb;
w=cb*inv（B）;
fori=1:
length（nobase）
sigma（i）=w*nB（:
i）-ncb（i）;
end
[maxs,ind]=max（sigma）;
ifmaxs<=0
minf=cb*xb;
vr=find（c~=0,1,'last'）;
forl=1:
vr
ele=find（baseVector==1,1）;
if（isempty（ele））
x（l）=0;
else
x（l）=xb（ele）;
end
end
bCon=0;
else
y=inv（B）*A（:
nobase（ind））;
ify<=0;
disp（'²»´æÔÚ×îÓŽâ'）;
x=NaN;
minf=NaN;
return;
else
minb=inf;
chagB=0;
forj=1:
length（y）
ify（j）>0
bz=xb（j）/y（j）;
ifbzminb=bz;
chagB=j;
end
end
end
tem=baseVector（chagB）;
baseVector（chagB）=nobase（ind）;
nobase（ind）=tem;
end
end
M=M+1;
if（M==10000）
disp（'ÕÒ²»µ½×îÓŽ⣡'）;
x=NaN;
minf=NaN;
return;
end
end
输入
>>A=[94100;310010;45001];
>>c=[-7-14000];
>>b=[360;300;200];
>>[x,mf]=simpleMthd（A,c,b,[3,4,5]）
输出结果
x=
2020
mf=
-476.0000
其他方法在实验报告中会有体现。
升级会员 