共轭梯度法
-
2023年3月17日发(作者：小羊过桥的故事)1
应用共轭梯度法求解方程组12
12
4
2
xx
xx
的根。初始值(0)(0,0)Tx
分析：将方程组
Axb
转化为优化问题中的极值问题然后应用共轭梯度法进行求解。
令
()RxAxb
，只需求得x，使得()Rx取得最小值。则有：
min()Rx2min()RxminTRR
()()
()()
()
()
()
()2
TT
TT
TTTTTT
TTTTTTT
TTTTT
TTTT
RRAxbAxb
XAbAxb
XAAXXAbbAxbb
XAAXXAbbAxbb
XAAXbAXbAxbb
XAAXbAXbb
这是一个数
与
1
()
2
TTfxxGxbxc
比较
则有：2,()22TTTGAAgxAAAb
回到题目问题中，则对应有1
2
412
4012
,,
04444
x
Ggb
x
也就是应用共轭梯度法求点使22
1212
()2212420fxxxxx取得最小值。
程序清单：
#include
#include
doublea,b,s[2];/*全局变量*/
doubleE=1e-6;
doubleF(doublex1,doublex2)
{
doubley;
y=2*x1*x1+2*x2*x2-12*x1-4*x2+20;
return(y);
}
voidqujian(doublex1,doublex2)/*用前进后退法求a的探索区间*/
{
doublea0=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)
{
2
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;} } } } doubleMIN(doublex1,doublex2)/*二次插值法求a的最小值*/ { doublex01,x02,x03,x0,xmin; doublef1,f2,f3,f; doubleX1,X2,X3,X4,X5,X6,X7,X8; doubleh=(b-a)/2; doublek1,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; 3 } elseif(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); } voidmain() { doubleg[3],m1,m2,b; intn=2,k=0; doublex1=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("最优解：nX1=%fnX2=%fn",x1,x2); } } 运行结果： 最优解： X1=3.000000 X2=1.000000 4