b=t2;t2=t1;f2=f1;
t1=*(b-a);f1=f(t1);
else
a=t1;t1=t2;f1=f2;
t2=a+*(b-a);f2=f(t2);
end
t=*(b+a);
k=k+1;
f0=f(t);
end
fprintf(1,'
迭代次数k=%\n',k)
fprintf(1,'
迭代区间一左端a=%\n',a)
fprintf(1,'
试点1坐标值t1=%\n',t1)
fprintf(1,'
函数值f1=%\n',f(t1))
fprintf(1,'
迭代区间一右端b=%\n',b)
fprintf(1,'
试点2坐标值t2=%\n',t2)
fprintf(1,'
函数值f2=%\n',f(t2))
fprintf(1,'
区间中点t=%\n',t)
fprintf(1,'
函数值fO=%\n',f(t))
运行结果如下:
迭代次数k=13
迭代区间一左端a=
试点1坐标值t1=
函数值f1=
迭代区间一右端b=
试点2坐标值t2=
函数值f2=
区间中点t=
函数值f0=
由黄金分割法在初始区间[0,3]求得的极小值点为t=,极小值为。
2
的极小点。
3.用牛顿法、阻尼牛顿法及变尺度法求函数fx1,x2x124x12x2
(1)在用牛顿法在MATLAB的M文件编辑器中编写的M文件,如下:
function[x,fx,k]=niudunfa(x0)symsx1x2
f=(x1-2)M+(x1-2*x2)A2;
fx=0;
v=[x1,x2];
df=jacobian(f,v);
df=df.';
G=jacobian(df,v);
epson=1e-12;
g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
G1=subs(G,{x1,x2},{x0(1,1),x0(2,1)});
k=0;
p=-G1\g1;
x0=x0+p;
while(norm(g1)>epson)
p=-G1\g1;
x0=x0+p;
g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
G1=subs(G,{x1,x2},{x0(1,1),x0(2,1)});
k=k+1;
end
x=x0;
fx=subs(f,{x1,x2},{x(1,1),x(2,1)});
运行结果如下:
>>[x,fx,k]=niudunfa([1;1])
x=
fx=
k=23
(2)用阻尼牛顿法在MATLA啲M文件编辑器中编写的M文件,如下:
function[x,fx,k]=zuniniudunfa(x0)%阻尼牛顿法
symsx1x2
f=(x1-2)A4+(x1-2*x2)A2;
fx=0;
v=[x1,x2];
df=jacobian(f,v);
df=df.';
G=jacobian(df,v);epson=1e-12;%停机原则
g1=subs(df,{x1,x2},{x0(1,1),x0(2,1)});
G1=subs(G,{x1,x2},{xO(1,1),xO(2,1)});
k=0;%迭代次数
p=-G1\g1;
aO=-p'*g1/(p'*G1*p);
xO=xO+aO*p;
while(norm(aO*p)>epson)
p=-G1\g1;
a0=-p'*g1/(p'*G1*p);
xO=xO+aO*p;
g仁subs(df,{x1,x2},{x0(1,1),x0(2,1)});
G1=subs(G,{x1,x2},{xO(1,1),xO(2,1)});
k=k+1;
end
x=x0;
fx=subs(f,{x1,x2},{x0(1,1),x0(2,1)});
运行结果如下:
>>[x,fx,k]=zuniniudunfa([1;1])
x=
fx=
k=23
(3)用变尺度法在MATLAB的M文件编辑器中编写的M文件,如下:
4.用共轭梯度法求函数fX「X2X;X;X/22为的极小点
22
(1)用共轭梯度法在MATLAB勺M文件编辑器中编写的M文件,如下:
function[y,x,k]=CG(A,b,c,xO)
%共轭梯度法解minf(x)=*X'*A*X+b'x+c
eps=1e-6;%迭代停机原则
%fx=*x0'.*A.*x0+b'.*x0+c;
r0=A*x0+b;
ifnorm(r0)<=eps
x=x0;
y=*x'*A*x+b'*x+c;
k=0;
end
p0=-r0;
a=-r0'*p0/(p0'*A*p0);
x1=x0+a*p0;
r1=A*x1+b;
k=0;
whilenorm(r1)>eps
beta=(r1'*r1)/(r0'*r0);p1=-r1+beta*p0;alpha=-(r1'*p1)/(p1'*A*p1);
x1=x1+alpha*p1;r2=A*x1+b;
p0=p1;
r0=r1;
r1=r2;
k=k+1;endx=x1;y=*x'*A*x+b'*x+c;运行结果如下:
[y,x,k]=CG([3-1;-11],[-2;0],0,[2;1])y=-1x=
M文件,如下:
k=1
(2)用变尺度法在MATLAB的M文件编辑器中编写的function[x,fx,k]=bianchidufa(A,b,c,x0)%用变尺度法求fx=*x'*A*x+b'*x+c;
epson=1e-12;g0=A*x0+b;
G0=A;H0=eye
(2);
k=0;d0=-H0*g0;a0=-d0'*g0/(d0'*G0*d0);s0=a0*d0;%x(k+1)-x(k);y0=A*a0*d0;%g(k+1)-g(k);x1=x0+a0*d0;
while(norm(s0)>=epson)
switchk
case{10}
x0=x1;g0=A*x0+b;H0=eye
(2);k=0;
d0=-H0*g0;a0=-d0'*g0/(d0'*G0*d0);
s0=a0*d0;x1=x0+a0*d0;
break
otherwise
g1=A*x1+b;
y0=A*a0*d0;
s0=a0*d0;
%H1=H0+s0*s0'/(s0'*y0)-H0*y0*y0'*H0/(y0'*H0*y0);
H1=H0+((1+y0'*H0*y0/(s0'*y0))*s0*s0'-H0*y0*s0'-s0*y0'*H0)/(s0'*y0);
k=k+1;
d1=-H1*g1;a1=-d1'*g1/(d1'*G0*d1);a0=a1;
d0=d1;
H0=H1;
s0=a0*d0;x1=x1+a0*d0;
break
end
end
x=x1;
fx=*x1'*A*x1+b'*x1+c;运行结果如下:
-1;-11],[-2;0],0,[2;1])
》[x,fx,k]=bianchidufa([3
H1=fx=-1x=
fx=-1
k=1
故函数极小点是点(1,1)
5.用鲍威尔法求函数fx1,x2x122x224x12x1x2的极小点。
用鲍威尔法在MATLAB的M文件编辑器中编写的M文件,如下:
function[x,fx,k]=bowell(A,b,c,x0)%鲍威尔法
d01=[1;0];
d02=[0;1];x02=[0;0];
esp=1e-12;%停机原则
k=0;%迭代次数
whilenorm(x0-x02)>=esp
k=k+1;
g01=A*x0+b;a01=-d01'*g01/(d01'*A*d01);
x01=x0+a01*d01;
g02=A*x01+b;
a02=-d02'*g02/(d02'*A*d02);
x02=x01+a02*d02;
d10=x02-x0;
g10=A*x02+b;
a10=-d10'*g10/(d10'*A*d10);
x10=x0+a01*d01;
d01=d02;
d02=d10;
x0=x10;
end
x=x0;
fx=*x'*A*x+b'*x+c;
运行结果如下:
[x,fx,k]=bowell([2-2;-24],[-4;0],0,[2;1])fx=
-8
x=
4
2
fx=
-8
k=
3
6.用单纯形法求线性规划问题
minf(x)1.1x12.2x23.3x34.4x4
x1x2x34
s.t.x12x2
2.5x33x45
xj0(j1,2,3,4)
用单纯形法在MATLAB的M文件编辑器中编写的M文件,如下:
%单纯形法matlab程序-danchunxingfa
%求解标准型线性规划:
maxc*x;.A*x=b;x>=0
列是资源向量b
0
%本函数中的A是单纯初始表,包括:
最后一行是初始的检验数,最后
%N是初始的基变量的下标
%输出变量sol是最优解,其中松弛变量(或剩余变量)可能不为
%输出变量val是最优目标值,kk是迭代次数function[sol,val,kk]=danchunxingfa(A,N)[mA,nA]=size(A);kk=0;%迭代次数flag=1;
whileflag
kk=kk+1;
ifA(mA,:
)<=0%已找到最优解
flag=0;
sol=zeros(1,nA-1);
fori=1:
mA-1
sol(N(i))=A(i,nA);
end
val=-A(mA,nA);
else
fori=1:
nA-1
ifA(mA,i)>0&A(1:
mA-1,i)<=0%问题有无界解disp('haveinfinitesolution!
');flag=0;
break;
end
end
ifflag%还不是最优表,进行转轴运算
temp=0;
fori=1:
nA-1
ifA(mA,i)>temp
temp=A(mA,i);
inb=i;%进基变量的下标
end
end
sita=zeros(1,mA-1);
fori=1:
mA-1
ifA(i,inb)>0
sita(i)=A(i,nA)/A(i,inb);
end
end
temp=inf;
fori=1:
mA-1
ifsita(i)>0&sita(i)temp=sita(i);
outb=i;%出基变量下标
end
end
%以下更新N
fori=1:
mA-1
ifi==outb
N(i)=inb;
end
end
%以下进行转轴运算
A(outb,:
)=A(outb,:
)/A(outb,inb);
ifi~=outb
A(i,:
)=A(i,:
)-A(outb,:
)*A(i,inb);
end
end
end
end
end;
令g(x)1.1x2.2x23.3x34.4x4,则求minf(x)1.1x12.2x23.3x34.4x4
就变成求maxg(x),即minf(x)maxg(x).
运行结果如下:
>>A=[11104;
1235;
0];N=[3;4];[sol,val,kk]=danchunxingfa(A,N)
sol=
00
val=
kk=
2
所以,
7.3333
经两次转轴运算,得到的最优解为x1x20,x34.0000,x41.667,minf(x)
7.求解线性规划问题
minz7x112x2
9x14x2x3360
4x15x2x4200
s.t.
3x110x2x5300xj0(j1,2,3,4,5)
用单纯形法在MATLAB的M文件编辑器中编写的M文件,如下:
%单纯形法matlab程序-danchunxingfa
%求解标准型线性规划:
maxc*x;.A*x=b;x>=0
%本函数中的A是单纯初始表,包括:
最后一行是初始的检验数,最后一列是资源向量b%N是初始的基变量的下标
%输出变量sol是最优解,其中松弛变量(或剩余变量)可能不为0
%输出变量val是最优目标值,kk是迭代次数function[sol,val,kk]=danchunxingfa(A,N)
[mA,nA]=size(A);
kk=0;%迭代次数
flag=1;
whileflag
kk=kk+1;
ifA(mA,:
)<=0%已找到最优解flag=0;
sol=zeros(1,nA-1);
fori=1:
mA-1sol(N(i))=A(i,nA);endval=-A(mA,nA);else
fori=1:
nA-1
ifA(mA,i)>0&A(1:
mA-1,i)<=0%问题有无界解disp('haveinfinitesolution!
');flag=0;
break;
end
end
ifflag%还不是最优表,进行转轴运算temp=0;
fori=1:
nA-1
ifA(mA,i)>temptemp=A(mA,i);inb=i;%进基变量的下标end
endsita=zeros(1,mA-1);
fori=1:
mA-1ifA(i,inb)>0sita(i)=A(i,nA)/A(i,inb);
end
endtemp=inf;
fori=1:
mA-1
ifsita(i)>0&sita(i)outb=i;%出基变量下标
end
end
%以下更新N
fori=1:
mA-1ifi==outbN(i)=inb;
end
end
%以下进行转轴运算A(outb,:
)=A(outb,:
)/A(outb,inb);
fori=1:
mA
ifi~=outb
A(i,:
)=A(i,:
)-A(outb,:
)*A(i,inb);
end
end
end
end
end;
运行结果如下:
令Y7xi12x2,则求minz7xi12x2可转变为求maXY,即minz-maXY.
>>A=[94100360;
45010200;
310001300;
7120000];N=[3;4;5];[sol,val,kk]=danchunxingfa(A,N)
sol=
20248400
val=
420
kk=
3所以,经3次转轴运算,得到的最优解为
x120,x224,x384,x4x50,minz420.