共轭梯度法求极小值

应用共轭梯度法求解方程组121242

x x x x +=⎧⎨-=⎩的根。初始值(0)(0,0)T x = 分析：将方程组Ax b =转化为优化问题中的极值问题然后应用共轭梯度法进行求解。 令()R x Ax b =-，只需求得x ，使得()R x 取得最小值。则有： min ()R x ⇔2min ()R x min T R R ⇔

()()

()()

()()()()2T T T T T T T T T T T T T T T T T T T T T T T T T T R R Ax b Ax b X A b Ax b X A AX X A b b Ax b b

X A AX X A b b Ax b b

X A A X b AX b Ax b b

X A A X b AX b b

=--=--=--+=--+=--+=-+这是一个数 与1()2

T T f x x Gx b x c =++比较 则有：2,()22T T T G A A g x A A A b ==- 回到题目问题中，则对应有124124012,,04444x G g b x --⎛⎫⎛⎫⎛⎫=== ⎪ ⎪ ⎪--⎝⎭⎝⎭⎝⎭

也就是应用共轭梯度法求点使221212()2212420f x x x x x =+--+取得最小值。

程序清单：

#include

#include

double a,b,s[2]; /*全局变量*/

double E=1e-6;

double F(double x1,double x2)

{

double y;

y=2*x1*x1+2*x2*x2-12*x1-4*x2+20;

return (y);

}

void qujian(double x1,double x2) /*用前进后退法求a 的探索区间*/ {

double a0=0,h=1,a1,a2,f1,f2,X1,X2,X3,X4;

a1=a0;a2=a0+h;

X1=x1+a1*s[0];

X2=x2+a1*s[1];

f1=F(X1,X2);

X3=x1+a2*s[0];

X4=x2+a2*s[1];

f2=F(X3,X4);

if (f1>f2)

{

while(1)

{

h=h*2;

a2=a2+h;

f1=f2;

X3=x1+a2*s[0];

X4=x2+a2*s[1];

f2=F(X3,X4);

if (f1>f2) {a1=a2-h;}

else{a=a1;b=a2;break;}

}

}

else

{

h=-h/4;

while(1)

{

a1=a1+h;

f2=f1;

X1=x1+a1*s[0];

X2=x2+a1*s[1];

f1=F(X1,X2);

if(f1

else{a=a1;b=a2;break;}

}

}

}

double MIN(double x1,double x2) /*二次插值法求a的最小值*/ {

double x01,x02,x03,x0,xmin;

double f1,f2,f3,f;

double X1,X2,X3,X4,X5,X6,X7,X8;

double h=(b-a)/2;

double k1,k2;

x01=a;x02=a+h;x03=b;

X1=x1+x01*s[0];

X2=x2+x01*s[1];

X3=x1+x02*s[0];

X4=x2+x02*s[1];

X5=x1+x03*s[0];

X6=x2+x03*s[1];

f1=F(X1,X2);

f2=F(X3,X4);

f3=F(X5,X6);

while (1)

{

k1=(f3-f1)/(x03-x01);

k2=((f2-f1)/(x02-x01)-k1)/(x02-x03);

x0=(x01+x03-k1/k2)/2;

X7=x1+x0*s[0];

X8=x2+x0*s[1];

f=F(X7,X8);

if(fabs(x0-x02)<=E)

{

xmin=x0;break;

}

else if(x02

{

if(f2

{x03=x0;f3=f;}

else

{x01=x02;f1=f2;x02=x0;f2=f;}

}

else

{

if(f2

{x01=x0;f1=f;}

else

{x03=x02;f3=f2;x02=x0;f2=f;}

}

}

return(xmin);

}

void main()

{

double g[3],m1,m2,b;

int n=2,k=0;

double x1=0,x2=0; /*起始点*/

g[1]=4*x1-12;

g[2]=4*x2-4;

m1=g[1]*g[1]+g[2]*g[2];

s[0]=-g[1];

s[1]=-g[2];

while(1)

{

qujian(x1,x2);

a=MIN(x1,x2); /*求a的探索区间*/

x1=x1+a*s[0];

x2=x2+a*s[1];

g[1]=4*x1-12;

m2=g[1]*g[1]+g[2]*g[2];

if(sqrt(g[1]*g[1]+g[2]*g[2])<=E)

break;

if(k+1==n)

{g[1]=4*x1-12;g[2]=4*x2-4;s[0]=-g[1];s[1]=-g[2];}

else

{b=m2/m1;s[0]=-g[1]+b*s[0];s[1]=-g[2]+b*s[1];k=k+1;}

printf("最优解：\n X1=%f\n X2=%f\n",x1,x2);

}

}

运行结果：

最优解：

X1=3.000000

X2=1.000000