共轭梯度法,黄金分割,外推法C语言代码

#include
#include
#define lamda 0.618
//************************
float a,b,c,d,e,g;//目标函数的系数
float x1,x2,arf;//目标函数的变量
float a0[2],b0[2],c0[2];//一维搜索区间的变量定义
float xi1,xi2,xi3,xi4;
float min_gold_x[2];
float directon[2];//搜索方向
float powell_minx[2],xg_min[2];
float _g[2];


//***********************************
void get_num()
{
printf("输入目标函数系数\n");
scanf("%f %f %f %f %f %f",&a,&b,&c,&d,&g,&e);
printf("a=%f b=%f c=%f d=%f e=%f g=%f\n",a,b,c,d,e,g);
printf("\n目标函数为：\n");
printf("y=%f*x1*x1+%f*x2*x2+%f*x1+%f*x2+%f*x1*x2+%f\n",a,b,c,d,g,e);
}

//功能函数***************************
float fun(float x1,float x2)
{
float y;
y=a*x1*x1+b*x2*x2+c*x1+d*x2+g*x1*x2+e;
return y;
}
//*********计算步长***
float change(float h0,float d[2])
{
float m,n;
n=(d[0]*d[0]+d[1]*d[1])/(h0*h0);
m=sqrt(n);
return m;
}

//***区间搜索********************
void search(float x01,float x02,float h0,float d[2])
{
float y0,y1,y2,n,di[2],x1[2],x2[2],x0[2];
int j,k=0;
for(j=0;j<2;j++)//读取搜索方向
{
di[j]=d[j];
}
x0[0]=x01,x0[1]=x02;
n=change(h0,di);
di[0]=di[0]/n;
di[1]=di[1]/n;

x1[0]=x0[0]+di[0];
x1[1]=x0[1]+di[1];

y0=fun(x0[0],x0[1]);
y1=fun(x1[0],x1[1]);
if(y1<=y0)
{
for(j=0;j<2;j++)
{
di[j]=2*di[j];
x2[j]=x1[j]+di[j];
}
}
if(y1>y0)
{
for(j=0;j<2;j++)
{
di[j]=(-2)*di[j];
x2[j]=x1[j]+di[j];
}
}
while(1)
{
y2=fun(x2[0],x2[1]);
if(y2<=y1)
{
y0=y1;y1=y2;
for(j=0;j<2;j++)
{
x0[j]=x1[j];
x1[j]=x2[j];
}
for(j=0;j<2;j++)
{
di[j]=2*di[j];
x2[j]=x1[j]+di[j];
}
}
if(y2>y1)
{
for(j=0;j<2;j++)
{
a0[j]=x0[j];
b0[j]=x1[j];
c0[j]=x2[j];
}
break;
}
}
}


//*********************黄金分割********************
void gold_d(float a[2],float b[2],float dlt)
{
float alf1[2],alf2[2],y1,y2;
float at[2],bt[2],temp1[2],temp2[2];
float m,n,i,k;
int j;
/*dlt=dlt*dlt;*/
for(j=0;j<2;j++)
{
at[j]=a[j];
bt[j]=b[j];
}
for(j=0;j<2;j++)
{
alf1[j]=bt[j]-lamda*(bt[j]-at[j]);
alf2[j]=at[j]+lamda*(bt[j]-at[j]);
}
y1=fun(alf1[0],alf1[1]);
y2=fun(alf2[0],alf2[1]);
while(1)
{
if(y1>=y2)
{
for(j=0;j<2;j++)
{
at[j]=alf1[j];
alf1[j]=alf2[j];
}
y1=y2;
for(j=0;j<2;j++)
{
alf2[j]=at[j]+lamda*(bt[j]-at[j]);
}
y2=fun(alf2[0],alf2[1]);
}
else
{
for(j=0;j<2;j++)
{
bt[j]=alf2[j];
alf2[j]=alf1[j];
}
y2=y1;
for(j=0;j<2;j++)
{
alf1[j]=bt[j]-lamda*(bt[j]-at[j]);
}
y1=fun(alf1[0],alf1[1]);
}

for(j=0;j<2;j++)
{
temp1[j]=bt[j]-at[j];
}
n=(y2-y1)/y2;n=fabs(n);
m=temp1[0]*temp1[0]+temp1[1]*temp1[1];

k=m/(bt[0]*bt[0]+bt[1]*bt[1]);
k=sqrt(k);m=sqrt(m);
i=y2-y1;i=fabs(i);

if((k{
for(j=0;j<2;j++)
{
min_gold_x[j]=(at[j]+bt[j])/2;
}
break;
}
}
}

/*******梯度********/
void g0(float x[2])
{
_g[0]=2*a*x[0]+c+g*x[1];
_g[1]=2*b*x[1]+d+g*x[0];
}
//*********************共轭梯度函数***********************************
void getd(float x0[2],float h,float dlt)
{

int i,j,k=0,q=0;
float t1,t2,t3,blt;
float dd[2],xt[2],xt1[2];
t1=100;
for(i=0;i<2;i++)
{
xt[i]=x0[i];
}


while(t1>=dlt&&q<5)
{
g0(xt);
for(i=0;i<2;i++)
{
dd[i]=-_g[i]; printf("\ng[0]=%f g[1]=%f\n",_g[0],_g[1]);
}
for(j=0;j<2;j++)
{

q++;
printf("\n第%d次\n",q);
search(xt[0],xt[1],h,dd);
gold_d(a0,c0,0.00001);
for(i=0;i<2;i++)
{
xt1[i]=min_gold_x[i]; printf("\nmin_gold_x[0]=%f min_gold_x[1]=%f\n",min_gold_x[0],min_gold_x[1]);
}
g0(xt1); printf("\ng[0]=%f g[1]=%f\n",_g[0],_g[1]);
t1=_g[0]*_g[0]+_g[1]*_g[1]; printf("\nt1=%f\n",t1);
t2=_g[0]*_g[0]+_g[1]*_g[1]; printf("\n|gk+1|2=%f\n",t2);
t3=dd[0]*dd[0]+dd[1]*dd[1]; printf("\n|gk|=%f\n",t3);
blt=t2/t3; printf("\nblt=%f\n",blt);
for(i=0;i<2;i++)
{
dd[i]=-_g[i]+blt*dd[i]; printf("\ndd[%d]=%f\n",i,dd[i]);
}
for(i=0;i<2;i++)
{
xt[i]=xt1[i];
}
}
for(i=0;i<2;i++)
{
xg_min[i]=xt1[i];
}
for(i=0;i<2;i++)
{
xt[i]=xt1[i];
}
}

}

//***********************主函数**************************
void main()
{
float x1[2];
float hh,z,y;
get_num();
printf("输入初始点坐标：\n");
scanf("%f %f",&x1[0],&x1[1]);
printf("初始点坐标为%f %f\n",x1[0],x1[1]);
printf("输入初始步长：\n");
scanf("%f",&hh);
printf("初始步长为：%f\n",hh);
printf("输入精度：\n");
scanf("%f",&z);
printf("要求精度为：%f\n",z);
getd(x1,hh,z);
y=fun(xg_min[0],xg_min[1]);
printf("\nxg_min[1]=%f xg_min[2]=%f\n",xg_min[0],xg_min[1]);
printf("y=%f",y);
}