共轭梯度算法

#include<stdio.h>
#include<math.h>
double x0[2],x1[2],z1,z2,z3,df[2],d1[2],d4[2],d0[2],a,c,c1,c2,c3,b,gg0,gg1,dd;/*x为坐标点，z为函数值，df导函数，d方向，ab区间，c步长，gg梯度的平方*/
int i,k,n=2;
double func(double x[],double d[],double c)
{
double z;
z=60-10*(x[0]+d[0]*c)-4*(x[1]+d[1]*c)+(x[0]+d[0]*c)*(x[0]+d[0]*c)+(x[1]+d[1]*c)*(x[1]+d[1]*c)-(x [0]+d[0]*c)*(x[1]+d[1]*c);/*原函数*/
return z;
}
void dfunc(double x[])/*导函数*/
{
df[0]=-10+2*x[0]-x[1];
df[1]=-4+2*x[1]-x[0];
}
void waitui(double x[],double d[],double c)/*外推*/
{
double h=0.03; /*起始步长*/
c1=0; /*起始点*/
z1=func(x,d,c1);
c2=c1+h;
z2=func(x,d,c2);
if(z2>z1) /*开始的下一个点大于初始点就反响收索，且交换开始两点的值*/
{
h=-h; /*反响*/
c3=c1;z3=z1; /*相互交换*/
c1=c2;z1=z2;
c2=c3;z2=z3;
}
c3=c2+h;z3=func(x,d,c3);/*求第三个点*/
while(z3<z2) /*迭代收索直到出现凹图形*/ {
h=h*2; /*不满足加大步伐*/
c1=c2;z1=z2; /*为下次迭代准备*/
c2=c3;z2=z3;
c3=c2+h;z3=func(x,d,c3);
}
a=c1; /*输出边界啊，（a，b）*/
b=c3;
}
double golden(double x[],double d[],double c) /*黄金分割*/ {
double j=0.618;
waitui(x,d,c); /*确定区间*/
c1=b-j*(b-a);
c2=a+j*(b-a);
z1=func(x,d,c1);
z2=func(x,d,c2);
if(z1>=z2) /*迭代，求下一个点以及函数值*/
{
a=c1;
c1=c2;z1=z2;
c2=a+j*(b-a);
z2=func(x,d,c2);
}
else
{
b=c2;
c2=c1;z2=z1;
c1=b-j*(b-a);
z1=func(x,d,c1);
}
while(abs(z2-z1)>=0.000001&&abs(c2-c1)>=0.000001||abs((z2-z1)/z2)>=0.000001&&abs((c 2-c1)/c2)>=0.000001) /*终止判断*/
{
if(z1>=z2) /*迭代，求下一个点以及函数值*/
{
a=c1;
c1=c2;z1=z2;
c2=a+j*(b-a);
z2=func(x,d,c2);
}
else
{
b=c2;
c2=c1;z2=z1;
c1=b-j*(b-a);
z1=func(x,d,c1);
}
}
c1=0.5*(c1+c2);
return c1;
}
void main()
{ k=1; /*迭代次数赋值*/ d4[0]=0;d4[1]=0; /*求函数值所用*/
for(i=0;i<n;i++)x0[i]=0;
dfunc(x0); /*初始点的梯度*/
for(i=0;i<n;i++)d0[i]=-df[i]; /*求负方向*/
while(k<200)
{
printf("\n\n第k=%3d迭代\n",k++);
printf("初始点为\n");
for(i=0;i<n;i++)
printf("x0(%d)=%12.6f\n\n",i,x0[i]);
c=golden(x0,d0,c); /*沿d0方向求最优步长*/
gg0=0; /*刷新*/
for(i=0;i<n;i++)gg0+=d0[i]*d0[i]; /*求平方*/
for(i=0;i<n;i++)x1[i]=x0[i]+d0[i]*c; /*求终点的坐标*/
printf("终点点为\n");
for(i=0;i<2;i++)
printf("x1(%d)=%12.6f\n\n",i,x1[i]);
dfunc(x1); /*终点的梯度*/
for(i=0;i<n;i++)d1[i]=df[i];
gg1=0;
for(i=0;i<n;i++)gg1+=d1[i]*d1[i];
z1=func(x1,d4,0); /*终点的函数值*/
if(sqrt(gg1)<=0.0001) /*结束判断*/
break;
dd=gg1/gg0;
for(i=0;i<n;i++)d1[i]=-d1[i]+dd*d0[i]; /*下一方向*/
for(i=0;i<n;i++)d0[i]=d1[i]; /*赋初值*/
for(i=0;i<n;i++)x0[i]=x1[i]; /*赋初始方向*/
if(k>200) /*达到迭代次数从新搜索*/
{ k=0;
dfunc(x0);
for(i=0;i<n;i++)d0[i]=-df[i];
}
}
printf("最终点为\n");
for(i=0;i<2;i++)
printf("x1(%d)=%12.6f\n\n",i,x1[i]);
printf("最终函数值为z1=%12.6f",z1);
}。