鲍威尔法C源程序.docx
《鲍威尔法C源程序.docx》由会员分享，可在线阅读，更多相关《鲍威尔法C源程序.docx（9页珍藏版）》请在冰豆网上搜索。
鲍威尔法C源程序
#include
#include
#include
#include
#include
#include
#definen12
#definett0.005
#definead0.0000001
//定义常量
//tt 初始迭代步长
//ad 收敛精度
floatia;
floatfny（float*x）
{
floatf;
f=10*pow（（x[0]+x[1]-5）,2）+pow（（x[0]-x[1]）,2）;//目标函数
return（f）;
}
float*iterate（float*x,floata,float*s）
{
float*x1;
x1=（float*）malloc（n1*sizeof（float））;
for（inti=0;ix1[i]=x[i]+a*s[i];
return（x1）;
}
floatfunc（float*x,floata,float*s）
{
float*x1;
x1=iterate（x,a,s）;
floatf=fny（x1）;
return（f）;
}
voidfinding（floata[3],floatf[3],float*xk,float*s）
{
floatt=tt;
floata1,f1;
a[0]=0;f[0]=func（xk,a[0],s）;
for（inti=0;;i++）
{
a[1]=a[0]+t;f[1]=func（xk,a[1],s）;
if（f[1]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;
}
//secondinsert
floatlagrange（float*xk,float*ft,float*s）
{
floata[3],f[3];
floatb,c,d,aa;
finding（a,f,xk,s）;
for（inti=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（*ft>f[1]）aa=a[1];
break;
}
if（aa>a[1]）
{
if（*ft>f[1]）
{
a[2]=aa;f[2]=*ft;
}
elseif（*ft{
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（*ft>f[1]）
{
a[0]=aa;f[0]=*ft;
}
elseif（*ft{
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（*ft>f[1]）
{
aa=a[1];*ft=f[1];
}
return（aa）;
}
float*powell（float*xk）
{
floath[n1][n1],s[n1]={0,0},ff[n1+1]={0,0,0};
floatf1,f3,aa;
floatdk[n1],*x0,xk1[n1];
intm=0,i,j;
for（i=0;i{
for（j=0;j{
h[i][j]=0;
if（j==i）
h[i][j]=1;
}
}
for（intk=0;;k++）
{
ff[0]=fny（xk）;
x0=xk;
for（i=0;i{
for（j=0;js[j]=h[i][j];
floataa=lagrange（xk,&ff[i+1],s）;
xk=iterate（xk,aa,s）;
}
for（i=0;i{
floata,b;
dk[i]=xk[i]-x0[i];
xk1[i]=2*xk[i]-x0[i];
}
floatmax=fabs（ff[1]-ff[0]）;
for（i=1;iif（fabs（ff[i+1]-ff[i]）>max）
{
max=fabs（ff[i+1]-ff[i]）;
m=i;
}
f3=fny（xk1）;
if（（f3{
aa=lagrange（xk,&f1,dk）;
xk=iterate（xk,aa,dk）;
for（i=m;ifor（j=0;jh[i][j]=h[i+1][j];
for（j=0;jh[n1-1][j]=dk[j];
}
else
{
if（ff[n1]>=f3）xk=xk1;
}
floatxq=0;
for（i=0;ixq+=pow（（xk[i]-x0[i]）,2）;
if（xq<=ad）
break;
}
return（xk）;
}
voidmain（）
{
floatxk[n1]={0,0};//取初始点
float*xx;
xx=（float*）malloc（n1*sizeof（float））;
xx=powell（xk）;
floatff=fny（xx）;
cout<<"优化的结果为：
"<printf（"\n\nTheOptimalDesignResultIs:
\n"）;
for（inti=0;iprintf（"\n\tx[%d]*=%f",i+1,xx[i]）;
printf（"\n\tf*=%f",ff）;
getch（）;
}
升级会员 