共轭梯度法求极小值

应用共轭梯度法求解方程组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 <stdio.h>
#include <math.h>
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<f2) {a2=a1-h;h=2*h;}
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<x0)
{
if(f2<f)
{x03=x0;f3=f;}
else
{x01=x02;f1=f2;x02=x0;f2=f;}
}
else
{
if(f2<f)
{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。