原文地址:. %powell·求方程 f=10*(x1+x2-5)^4+(x1-x2+x3)^2 +(x2+x3)^6 的最优解 f nction mypowell syms x1 x2 x3 a; f=10*(x1+x2-5)^4+(x1-x2+x3)^2 +(x2+x3)^6; error=10^(-3); d=eye(3); x0=[0 0 0] ; for k=1:1:10^6 maxlength=0;x00=x0;m=0; if k==1,s=d;end for i=1:3 x=x0+a*s(:,i); ff=s bs(f,{x1,x2,x3},{x(1),x(2),x(3)}); t=divide(ff,a); %调用了进退法分割区间 aa=onedemensionslsearch(ff,a,t); %调用了0.618法进行一维搜索 xx=x0+aa*s(:,i); fx0=s bs(f,{x1,x2,x3},{x0(1),x0(2),x0(3)}); fxx=s bs(f,{x1,x2,x3},{xx(1),xx(2),xx(3)}); length=fx0-fxx; if lengthmaxlength,maxlength=length;m=m+1;end x0=xx; end ss=x0-x00; reflectx=2*x0-x00; f1=s bs(f,{x1,x2,x3},{x00(1),x00(2),x00(3)}); f2=s bs(f,{x1,x2,x3},{x0(1),x0(2),x0(3)}); f3=s bs(f,{x1,x2,x3},{reflectx(1),reflectx(2),reflectx(3)}); if f3f1&&(f1+f3-2*f2)*(f1-f2-maxlength)^20.5*maxlength*(f1-f3)^2 x=x0+a*ss; ff=s bs(f,{x1,x2,x3},{x(1),x(2),x(3)}); t=divide(ff,a); aa=onedemensionslsearch(ff,a,t); x0=x0+aa*ss; for j=m:(3-1),s(:,j)=s(:,j+1);end s(:,3)=ss; else if f2f3, x0=reflectx;end end if norm(x00-x0)error,break;end k x0 end opx=x0; val=s bs(f,{x1,x2,x3},{opx(1),opx(2),opx(3)}); disp( 最优点: );opx disp( 最优化值: );val disp( 迭代次数: );k 附录: 0.618法 f nction o tp t=onedemensionslsearch(f,x,s,r) if nargin4,r=1e-6;end a=s(1);b=s(2); a1=a+0.382*(b-a);fa1=s bs(f,{x},{a1}); a2=a+0.618*(b-a);fa2=s bs(f,{x},{a2}); while abs((b-a)b)r && abs((fa2-fa1)fa2)r if fa1fa2 b=a2;a2=a1;fa2=fa1;a1=a+0.382*(b-a);fa1=s bs(f,{x},{a1}); else a=a1;a1=a2;fa1=fa2;a2=a+0.618*(b-a);fa2=s bs(f,{x},{a2}); end end op=(a+b)2; %fop=s bs(f,{x},{op}); o tp t=op; 进退法 %对任意一个一维函数函数进行区间分割,使其出现"高-低-高"的型式 f nction o tp t=divide(f,x,m,n) if nargin4,n=1e-6;end if nargin3,m=0;end step=n; t0=m;ft0=s bs(f,{x},{t0}); t1=t0+step;ft1=s bs(f,{x},{t1}); if ft0=ft1 t2=t1+step;ft2=s bs(f,{x},{t2}); while ft1ft2 t0=t1; %ft0=ft1; t1=t2;ft1=ft2; step=2*step;t2=t1+step;ft2=s bs(f,{x},{t2}); end else step=-step; t=t0;t0=t1;t1=t;ft=ft0; %ft0=ft1; ft1=ft; t2=t1+step;ft2=s bs(f,{x},{t2}); while ft1ft2 t0=t1; %ft0=ft1; t1=t2;ft1=ft2; step=2*step;t2=t1+step;ft2=s bs(f,{x},{t2}); end end o tp t=[t0,t2];