Powell法

Powell法的C语言程序程序使用说明：a[0]—初始区间的下界值；ad—一维迭代精度及收敛精度；tt—一维搜索的初始迭代步长；xk[n]—迭代初始点。程序如下：#includestdio.h#includestdlib.h#includemath.h#includealloc.h#includeconio.h#definen2#definett0.005#definead0.01doubleia;doublefny(double*x){doublef;f=pow(x[0],2)+2*pow(x[1],2)-4*x[0]-2*x[0]*x[1];returnf;}double*iterate(double*x,doublea,double*s){double*x1;inti;x1=(double*)malloc(n*sizeof(double));for(i=0;in;i++)x1[i]=x[i]+a*s[i];return(x1);}doublefunc(double*x,doublea,double*s){double*x1;doublef;x1=iterate(x,a,s);f=fny(x1);returnf;}voidfinding(doublea[3],doublef[3],double*xk,double*s){doublet=tt;doublea1,f1;inti;a[0]=0;f[0]=func(xk,a[0],s);for(i=0;;i++){a[1]=a[0]+t;f[1]=func(xk,a[1],s);if(f[1]f[0])break;if(fabs(f[1]-f[0])=ad){t=-t;a[0]=a[1];f[0]=f[1];}else{if(ia==1)return;t=t/2;ia=1;}}for(i=0;;i++){a[2]=a[1]+t;f[2]=func(xk,a[2],s);if(f[2]f[1])break;t=2*t;a[0]=a[1];f[0]=f[1];a[1]=a[2];f[1]=f[2];}if(a[0]a[2]){a1=a[0];f1=f[0];a[0]=a[2];f[0]=f[2];a[2]=a1;f[2]=f1;}return;}doublelagrange(double*xk,double*ft,double*s){inti;doublea[3],f[3];doubleb,c,d,aa;finding(a,f,xk,s);for(i=0;;i++){if(ia==1){aa=a[1];*ft=f[1];break;}d=(pow(a[0],2)-pow(a[2],2))*(a[0]-a[1])-(pow(a[0],2)-pow(a[1],2))*(a[0]-a[2]);if(fabs(d)==0)break;c=((f[0]-f[2])*(a[0]-a[1])-(f[0]-f[1])*(a[0]-a[2]))/d;if(fabs(c)==0)break;b=((f[0]-f[1])-c*(pow(a[0],2)-pow(a[1],2)))/(a[0]-a[1]);aa=-b/(2*c);*ft=func(xk,aa,s);if(fabs(aa-a[1])=ad){if(*ftf[1])aa=a[1];break;}if(aaa[1]){if(*ftf[1]){a[2]=aa;f[2]=*ft;}elseif(*ftf[1]){a[0]=a[1];a[1]=aa;f[0]=f[1];f[1]=*ft;}elseif(*ft==f[1]){a[2]=aa;a[0]=a[1];f[2]=*ft;f[0]=f[1];a[1]=(a[0]+a[2])/2;f[1]=func(xk,a[1],s);}}else{if(*ftf[1]){a[0]=aa;f[0]=*ft;}elseif(*ftf[1]){a[2]=a[1];a[1]=aa;f[2]=f[1];f[1]=*ft;}elseif(*ft==f[1]){a[0]=aa;a[2]=a[1];f[0]=*ft;f[2]=f[1];a[1]=(a[0]+a[2])/2;f[1]=func(xk,a[1],s);}}}if(*ftf[1]){aa=a[1];*ft=f[1];}return(aa);}double*powell(double*xk){doubleh[n][n],s[n]={0,0},ff[n+1]={0,0,0};doublef1,f3,max,aa=0.0,xq;double*dk,*x0,*xk1;intm=0,i,j,k,l;for(l=0;ln;l++){dk[l]=0;xk1[l]=0;}for(i=0;in;i++)for(j=0;jn;j++){h[i][j]=0;if(j==i)h[i][j]=1;}for(k=0;;k++){ff[0]=fny(xk);x0=xk;for(i=0;in;i++){for(j=0;jn;j++)s[j]=h[i][j];aa=lagrange(xk,&ff[i+1],s);xk=iterate(xk,aa,s);}for(i=0;in;i++){dk[i]=xk[i]-x0[i];xk1[i]=2*xk[i]-x0[i];}max=fabs(ff[1]-ff[0]);for(i=1;in;i++)if(fabs(ff[i+1]-ff[i])max){max=fabs(ff[i+1]-ff[i]);m=i;}f3=fny(xk1);if((f3ff[0])&&((ff[0]+f3-2*ff[n])*pow((ff[0]-ff[n]-max),2)0.5*max*pow((ff[0]-f3),2))){aa=lagrange(xk,&f1,dk);xk=iterate(xk,aa,dk);for(i=m;in-1;i++)for(j=0;jn;j++)h[i][j]=h[i+1][j];for(j=0;jn;j++)h[n-1][j]=dk[j];}else{if(ff[n]=f3)xk=xk1;}xq=0;for(i=0;in;i++)xq+=pow((xk[i]-x0[i]),2);if(xq=ad)break;}return(xk);}voidmain(){doublexk[n]={1,1};double*xx,ff;inti;clrscr();xx=(double*)malloc(n*sizeof(double));xx=powell(xk);ff=fny(xx);printf(\n\nTheOptimalDesignResultIs:\n);for(i=0;in;i++)printf(\n\tx[%d]=%f,i+1,xx[i]);printf(\n\tf=%f,ff);getch();}