鲍威尔法C++源程序

#include <iostream.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <malloc.h>
#include <conio.h>
#define n1 2
#define tt 0.005
#define ad 0.0000001
//定义常量
//tt初始迭代步长
//ad收敛精度
float ia;
float fny(float *x)
{
float f;
f=10*pow((x[0]+x[1]-5),2)+pow((x[0]-x[1]),2); //目标函数return(f);
}
float *iterate(float *x,float a, float *s)
{
float *x1;
x1=(float *)malloc(n1 * sizeof(float));
for (int i=0;i<n1;i++)
x1[i]=x[i]+a*s[i];
return(x1);
}
float func(float *x,float a,float *s)
{
float *x1;
x1=iterate(x,a,s);
float f=fny(x1);
return(f);
}
void finding(float a[3],float f[3],float *xk,float *s) {
float t=tt;
float a1,f1;
a[0]=0; f[0]=func(xk,a[0],s);
for (int 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;
}
//second insert
float lagrange(float *xk,float *ft,float *s)
{
float a[3],f[3];
float b,c,d,aa;
finding (a,f,xk,s);
for (int 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 (*ft>f[1]) aa=a[1];
break;
}
if (aa>a[1])
{
if (*ft>f[1])
{
a[2]=aa; f[2]=*ft;
}
else if (*ft<f[1])
{
a[0]=a[1]; a[1]=aa;
f[0]=f[1]; f[1]=*ft;
}
else if (*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;
}
else if (*ft<f[1])
{
a[2]=a[1]; a[1]=aa;
f[2]=f[1]; f[1]=*ft;
}
else if (*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)
{
float h[n1][n1],s[n1]={0,0},ff[n1+1]={0,0,0};
float f1,f3,aa;
float dk[n1],*x0,xk1[n1];
int m=0,i,j;
for (i=0;i<n1;i++)
{
for(j=0;j<n1;j++)
{
h[i][j]=0;
if (j==i)
h[i][j]=1;
}
}
for (int k=0;;k++)
{
ff[0]=fny(xk);
x0=xk;
for (i=0;i<n1;i++)
{
for (j=0;j<n1;j++)
s[j]=h[i][j];
float aa=lagrange(xk,&ff[i+1],s);
xk=iterate(xk,aa,s);
}
for (i=0;i<n1;i++)
{
float a,b;
dk[i]=xk[i]-x0[i];
xk1[i]=2*xk[i]-x0[i];
}
float max=fabs(ff[1]-ff[0]);
for (i=1;i<n1;i++)
if (fabs(ff[i+1]-ff[i])>max)
{
max=fabs(ff[i+1]-ff[i]);
m=i;
}
f3=fny(xk1);
if
((f3<ff[0])&&((ff[0]+f3-2*ff[2])*pow((ff[0]-ff[n1]-max),2)<0.5*max*pow((ff[0]-f3),2))) {
aa=lagrange(xk,&f1,dk);
xk=iterate(xk,aa,dk);
for (i=m;i<n1-1;i++)
for (j=0;j<n1;j++)
h[i][j]=h[i+1][j];
for (j=0;j<n1;j++)
h[n1-1][j]=dk[j];
}
else
{
if(ff[n1]>=f3) xk=xk1;
}
float xq=0;
for (i=0;i<n1;i++)
xq+=pow((xk[i]-x0[i]),2);
if (xq<=ad)
break;
}
return(xk);
}
void main()
{
float xk[n1]={0,0};//取初始点
float *xx;
xx=(float *)malloc(n1 *sizeof(float));
xx=powell(xk);
float ff=fny(xx);
cout<<"优化的结果为："<<endl;
printf("\n\nThe Optimal Design Result Is:\n");
for (int i=0;i<n1;i++)
printf("\n\t x[%d] *=%f",i+1,xx[i]);
printf("\n\t f*=%f",ff);
getch();
}。