机械优化设计MATLAB程序

t t t 机械优化设计作业

1.用二次插值法求函数ϕ( )= ( +1)( - 2)2 极小值，精度 e=0.01。

在 MA TLAB 的 M 文件编辑器中编写的 M 文件，如下：

f=inline('(t+1)*(t -2)^2','t')

a=0;b=3;epsilon=0.01;

t1=a;f1=f(t1);

t3=b;f3=f(t3);

t2=0.5*(t1+t3);f2=f(t2);

c1=(f3-f1)/(t3-t1);

c2=((f2-f1)/(t2-t1)-c1)/(t2-t3);

t4=0.5*(t1+t3-c1/c2);f4=f(t4);

k=0;

while(abs(t4-t2)>=epsilon)

if t2

if f2>f4

f1=f2;t1=t2;

t2=t4;f2=f4;

else

f3=f4;t3=t4;

end

else

if f2>f4

f3=f2;t3=t2;

t2=t4;f2=f4;

else

f1=f4;t2=t4;

end

end

c1=(f3-f1)/(t3-t1);

c2=((f2-f1)/(t2-t1)-c1)/(t2-t3);

t4=0.5*(t1+t3-c1/c2);f4=f(t4);

k=k+1;

end

%输出最优解

if f2>f4

t=t4;f=f(t4);

else

t=t2;f=f(t2);

end

fprintf(1,'迭代计算 k=%3.0f\n',k)

fprintf(1,'极小点坐标 t=%3.0f\n',t)

fprintf(1,'函数值 f=%3.4f\n',f)

3.用牛顿法、阻尼牛顿法及变尺度法求函数 的极小点。( ) ( ) ( )21121

22, xxxxxf -+-= 4 2 (1)在用牛顿法在 MATLAB 的 M 文件编辑器中编写的 M 文件，如下： function [x,fx,k]=niudunfa(x0)

syms x1 x2

f=(x1-2)^4+(x1-2*x2)^2;

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 =1.9999554476059523381489991377897

0.99997772380297616907449956889483

fx =0.0000000000000000039398907941382470301534502947647 k =23

(2)用阻尼牛顿法在 MA TLAB 的 M 文件编辑器中编写的 M 文件，如下： function [x,fx,k]=zuniniudunfa(x0)%阻尼牛顿法

syms x1 x2

f=(x1-2)^4+(x1-2*x2)^2;

fx=0;

v=[x1,x2];

df=jacobian(f,v);

df=df.';

G=jacobian(df,v);

epson=1e -12;%停机原则

k=0;

while norm(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;

end

x=x1;

y=0.5*x'*A*x+b'*x+c;

运行结果如下：

[y,x,k]=CG([3-1;-11],[-2;0],0,[2;1])

y=-1

x= 1.0000

1.0000

k=1

(2)用变尺度法在MATLAB的M文件编辑器中编写的M文件，如下：function[x,fx,k]=bianchidufa(A,b,c,x0)

%用变尺度法求fx=0.5*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)

switch k

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;