1、工程一班23号甘超第三次作业第三次作业P126 6-2一(1) 混合函数法所用到的牛顿法M文件functionx,minf=minNT(f,x0,var,eps)format long;if nargin=3 eps=1.0e-6;endtol=1;x0=transpose(x0);gradf=jacobian(f,var);jacf=jacobian(gradf,var);while toleps v=subs(gradf,var,x0); tol=norm(v); pv=subs(jacf,var,x0); p=-inv(pv)*transpose(v); p=double(p); x1=
2、x0+p; x0=x1;endx=x1;minf=subs(f,var,x);format short;所用到的混合函数M文件function x,minf=minMixFun6-2(f,g,x0,r0,c,var,eps)gx0=subs(g,var,x0);if gx0=0;else disp(!); x=NaN; minf=NaN; return;endif r0=1|c0 disp(01!); x=NaN; minf=NaN; return;endif nargin=7 eps=1.0e-6;endFE=0;for i=1:length(g) FE=FE+1/g(i);endFH=tr
3、anspose(h)*h;x1=transpose(x0);x2=inf;while 1 FF=r0*FE+FH/sqrt(r0); SumF=f+FF; x2,minf=minNT(SumF,transpose(x1),var); if norm(x2-x1) syms s t; f=s2+2*t2; 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.3333minf =0.3333二（2）内点罚函数法所插入的牛顿法M文件functionx,minf=minNT(f,x0,
4、var,eps)format long;if nargin=3 eps=1.0e-6;endtol=1;x0=transpose(x0);gradf=jacobian(f,var);jacf=jacobian(gradf,var);while toleps 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;endx=x1;minf=subs(f,var,x);format short;内点惩罚函数法M文件function x,
5、minf=minNF(f,x0,g,u,v,var,eps)format long;if nargin=6 eps=1.0e-6;endk=0;FE=0;for i=1:length(g) FE=FE+1/g(i);endx1=transpose(x0);x2=inf;while 1 FF=u*FE; SumF=f+FF; x2,minf=minNT(SumF,transpose(x1),var); Bx=subs(FE,var,x2); if u*Bxeps if norm(x2-x1)=eps x=x2; break; else u=u*v; x1=x2; endelse if norm(
6、x2-x1) clear syms t;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.0000minf = 50.0000单纯形法求解线性规划问题单纯形法M文件functionx,minf=simpleMthd(A,c,b,baseVector)sz=size(A);nVia=sz(2);n=sz(1);xx=1:nVia;nobase=zeros(1,1);m=1;for i=1:nVia; if (isempty(find(baseVector=xx(i),1) nobase(m)=i; m=m
7、+1; else ; endendbCon=1;M=0;while bCon nB=A(:,nobase); ncb=c(nobase); B=A(:,baseVector); cb=c(baseVector); xb=inv(B)*b; f=cb*xb; w=cb*inv(B); for i=1:length(nobase) sigma(i)=w*nB(:,i)-ncb(i); end maxs,ind=max(sigma); if maxs=0 minf=cb*xb; vr=find(c=0,1,last); for l=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); if y0 bz=xb(j)/y(j); if bz A=9 4 1 0 0;3 10 0 1 0; 4 5 0 0 1; c=-7 -14 0 0 0; b=360;300;200; x,mf=simpleMthd(A,c,b,3,4,5)输出结果x = 20 20mf = -476.0000其他方法在实验报告中会有体现。
