共轭梯度法源代码.docVIP

共轭梯度法程序源代码 #includestdio.h #includemath.h #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(f1f0) { h=alpha*h; t=t0; t0=t1; f0=f1; k++; } else { if(k==0) {h=-h;t=t1;} else { *a=tt1?t:t1; *b=tt1?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)eps) break; else { if(f1f2) { 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;in;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(mod1eps) { 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 { if(k+1==n) {