共轭梯度法源代码

共轭梯度法程序源代码

#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)); /*求梯度的长度*/

if(mod2<=eps) break;

else

{