共轭梯度法

共轭梯度法

1.算法思想：

共轭梯度法是利用目标函数梯度逐步产生共轭方向作为线搜索方向的方法，每次搜索方向都是在目标函数梯度的共轭方向，搜索步长通过一维极值算法确定。

2.算法步骤：

用共轭梯度法求无约束多维极值问题min (),n f x x R ∈的算法步骤如下：

（1） 给定初始点(0)x ，及精度0ε>； （2） 若

(0)()f x ε

∇≤，停止，极小值点为(0)x ，否则转步骤（3）；

（3） 取(0)(0)()p f x =-∇，且置0k =；

（4） 用一维搜索法求k t ，使得()()()()()0

()min

k k k k k t f x t p f x tp ≥+=+，令，(1)()()k k k k x x t p +=+，转步骤5； （5） 若

(1)()k f x ε

+∇≤，停止，极小值点为(1)k x +，否则转步骤（6）；

（6） 若1k n +=，令(0)()n x x =，转步骤（3），否则转步骤（7）； （7） 令(1)(1)()()k k k k p f x p λ++=-∇+，2(1)2

()

()

()

k k

k f x

f x λ+∇=∇，置1k k =+，转

步骤（4）。

3.算法源程序：

#include #include

#define N 10

#define eps pow(10,-6)

double f(double x[],double p[],double t)

{

double s;

s=pow(x[0]+t*p[0],2)+25*pow(x[1]+t*p[1],2); return s;

}

/*以下是进退法搜索区间源程序*/

void sb(double *a,double *b,double x[],double p[]) {

double t0,t1,t,h,alpha,f0,f1;

int k=0;

t0=2.5; /*初始值*/

h=1; /*初始步长*/

alpha=2; /*加步系数*/

f0=f(x,p,t0);

t1=t0+h;

f1=f(x,p,t1);

while(1)

{

if(f1

{

h=alpha*h; t=t0;

t0=t1; f0=f1;

k++;

}

else

{

if(k==0)

{h=-h;t=t1;}

else

{

*a=t

*b=t>t1?t:t1;

break;

}

}

t1=t0+h;

f1=f(x,p,t1);

}

}

/*以下是黄金分割法程序源代码*/

double hjfg(double x[],double p[])

{

double beta,t1,t2,t;

double f1,f2;

double a=0,b=0;

double *c,*d;

c=&a,d=&b;

sb(c,d,x,p);/*调用进退法搜索区间*/

printf("\nx1=%lf,x2=%lf,p1=%lf,p2=%lf",x[0],x[1],p[0],p[1]); printf("\n[a,b]=[%lf,%lf]",a,b);

beta=(sqrt(5)-1.0)/2;

t2=a+beta*(b-a); f2=f(x,p,t2);

t1=a+b-t2; f1=f(x,p,t1);

while(1)

{

if(fabs(t1-t2)

break;

else

{

if(f1

{

t=(t1+t2)/2;

b=t2; t2=t1;

f2=f1; t1=a+b-t2;

f1=f(x,p,t1);

}

else

{

a=t1; t1=t2;

f1=f2;

t2=a+beta*(b-a);

f2=f(x,p,t2);

}

}

}

t=(t1+t2)/2;

return t;

}

/*以下是共轭梯度法程序源代码*/

void gtd()

{

double x[N],g[N],p[N],t=0,f0,mod1=0,mod2=0,nanda=0; int i,k,n;

printf("请输入函数的元数值n=");

scanf("%d",&n);

printf("\n请输入初始值:\n");

for(i=0;i

scanf("%lf",&x[i]);

f0=f(x,g,t);

g[0]=2*x[0]; g[1]=50*x[1];

mod1=sqrt(pow(g[0],2)+pow(g[1],2));/*求梯度的长度*/

if(mod1>eps)

{

p[0]=-g[0]; p[1]=-g[1]; k=0;

while(1)

{

t=hjfg(x,p);/*调用黄金分割法求t的值*/

printf("\np1=%lf,p2=%lf,t=%lf",p[0],p[1],t);

x[0]=x[0]+t*p[0]; x[1]=x[1]+t*p[1];

g[0]=2*x[0]; g[1]=50*x[1];

/*printf("\nx1=%lf,x2=%lf,g1=%lf,g2=%lf",x[0],x[1],g [0],g[1]);*/

mod2=sqrt(pow(g[0],2)+pow(g[1],2)); /*求梯度的长度*/