广义最小残差法和共轭梯度法的区别

Posted 2023-04-09

答：共轭梯度法是以函数的梯度构造共轭方向的一种算法，具有共轭方向的性质。共轭梯度法具有超线性收敛速度。梯度法与共轭梯度法的区别是：
1)最速下降法(梯度法) ：搜索方向为目标函数负梯度方向，计算效率优于坐标轮换法。开始几步搜索下降快，但愈接近极值点下降愈慢。对初始点的选择要求不高，适合与其它方法结合使用。
2)共轭梯度法：第一步搜索沿负梯度方向，然后沿负梯度的共轭方向搜索。计算效率介于梯度法和牛顿法之间。对初始点没有特殊的要求，不需要计算二阶偏导数矩阵及其逆矩阵，计算量与梯度法相当。适用于各种规模的问题 参考技术A 【数值算法】系数矩阵非对称时，线性方程组如何求解？-稳定双共轭梯度法(Bicgstab)求解线性方程组

寒江雪
来自专栏有限元术
在前面的文章和中表明共轭梯度法是求解对称正定线性方程组的一种有效方法，当针对不同的系数矩阵采用不同的预处理方式时，其可以以较少的迭代次数获得较高精度的解。然而，该方法的一个缺点就是其只能适用于对称正定系数矩阵，当系数矩阵不再是对称正定时，此方法可能失效。以下举例：

上面矩阵A为非对称矩阵，采用共轭梯度法求解过程如下：

该方程组采用共轭梯度法迭代4862次依然未收敛。因此，对于该非对称方程，可以认为，共轭梯度法几乎是失效的。在实际工程中，有限元方法形成的刚度系数以对称正定居多，但是实际上也存在非对称的可能，例如，当材料本构采用摩尔-库伦本构时，其形成的刚度矩阵就有可能会是非对称的，此时如果是使用商业软件，应当在软件中选择非对称求解器。如果是自主编程且采用迭代法求解线性方程组，则需要找到一种适用于非对称矩阵的求解方法。

常见的非对称系数矩阵求解方法主要有：广义最小残差法(GMRES)，双共轭梯度法（Bicg）稳定双共轭梯度法(BiCGStab)，稳定混合双共轭梯度法（BiCGStab(l)），这些方法相对于常规的共轭梯度法在推导上均增加了一些难度，实际推导往往较为复杂。本文不展开推导，仅对稳定双共轭梯度法(BiCGStab)的伪代码作简要粘贴。BiCGStab法的具体计算过程如下：

具体代码：

program bicgstab_main implicit none integer,parameter::n=8 real(8)::a(n,n),b(n) real(8)::x(n),x0(n) integer::i,j integer,parameter::m=20 real(8)::c(m,m),d(m) real(8)::y(m),y0(m) a=reshape((/6,5,1,2,0,0,2,1, & 0,5,1,1,0,0,3,0,& 1,1,623,1,2,0,1001,2, & 2,1,1,7,1,2,1,1,& 0,0,2,31,6,0,2,1,& 0,0,0,2,0,4,1,0,& 2,3,1,1,23,1,5,213,& 123,0,12,1,1,0,1,3/),(/n,n/)) b=(/1,1,1,1,1,1,1,1/) write(*,*)"矩阵A" write(*,"(8f10.4)")(a(i,:),i=1,n) write(*,*)"向量b" write(*,"(f10.4)")b x0=100.0 x=0.0 call bicgstab(a,b,x,x0,n) end program bicgstab_main subroutine bicgstab(a,b,x,x0,n) implicit none integer,intent(in)::n real(kind=8),intent(in)::a(n,n),b(n),x0(n) real(kind=8),intent(out)::x(n) real(kind=8)::p0(n),r0(n) real(kind=8)::tol=1.0d-6 integer::nmax=1000 real(kind=8)::rj(n),pj(n) real(kind=8)::alphaj real(kind=8)::r0_bar(n) real(kind=8)::sj(n) real(kind=8)::xjp1(n),xj(n) real(kind=8)::wj real(kind=8)::rjp1(n) real(kind=8)::betaj integer::n_iter real(kind=8)::apj(n),asj(n) r0=b-matmul(a,x0) r0_bar=r0 if(abs(dot_product(r0,r0_bar))<tol) then write(*,*)"unvalid r0_bar,select an other!" return endif p0=r0 rj=r0 pj=p0 xj=x0 n_iter=0 do if(n_iter>nmax) then write(*,*)"exceed max iter!" exit endif alphaj=dot_product(rj,r0_bar)/dot_product(matmul(a,pj),r0_bar) apj=matmul(a,pj) sj=rj-alphaj*apj if(norm2(sj)<tol) then xjp1=xj+alphaj*pj x=xjp1 exit endif asj=matmul(a,sj) wj=dot_product(asj,sj)/(norm2(asj))**2 xjp1=xj+alphaj*pj+wj*sj rjp1=sj-wj*asj if(norm2(rjp1)<tol) then x=xjp1 exit endif betaj=alphaj*dot_product(rjp1,r0_bar)/(wj*dot_product(rj,r0_bar)) pj=rjp1+betaj*(pj-wj*apj) xj=xjp1 rj=rjp1 n_iter=n_iter+1 write(*,*)"the",n_iter,"iter!" enddo write(*,*)"the solution of equation:" write(*,"(es18.8)")x end subroutine bicgstab
依据上述过程编写程序，计算前述非对称矩阵线性方程组求解结果：

采用matlab求解该方程组的解：

通过对比可知11次迭代已经获得即为准确的结果。实际上，对于该方法也可以通过一定的预处理方式，使得其所需要的迭代次数更少。以上，就是稳定双共轭梯度法求解线性方程组的内容，感谢您的阅读！欢迎关注公众号 有限元术

请问有人知道共轭梯度法的FR,PRP,HS三个算法的Matlab程序吗?

如题...感激不尽!

%共轭梯度法 FR
% G为对称正定矩阵，X是初始点，e为精度
%a是精确线搜索步长
function [m2,a,d,X,g1,f1] = conjgrad(G,b,c,X,e)
n=length(G);
if n==2
format long e %rat
syms x1 x2
f=1/2*[x1,x2]*G*[x1;x2]+b'*[x1;x2]+c;
g=[diff(f,x1);diff(f,x2)];
g1=subs(subs(g,x1,X(1,1)),x2,X(2,1));
d=-g1;
a=-(d'*g1)/(d'*G*d);% a=-((X(:,1)'*G*d+b'*d)/(d'*G*d)); a=g1(:,1)'*g1(:,1)/(d(:,1)'*G*d(:,1));
X(:,2)=X(:,1)+a*d;
g1=[g1 subs(subs(g,x1,X(1,2)),x2,X(2,2))];
m1=norm(g1(:,1));m2=norm(g1(:,2));
k=(m2/m1)^2;
i=2;
while m2>=e
d(:,i)=-g1(:,i)+k(i-1)*d(:,i-1);
a(i)=-(d(:,i)'*g1(:,i))/(d(:,i)'*G*d(:,i));
%a1(i)=-((X(:,i)'*G*d(:,i)+b'*d(:,i))/(d(:,i)'*G*d(:,i)));a(i)=g1(:,i)'*g1(:,i)/(d(:,i)'*G*d(:,i));
X(:,i+1)=X(:,i)+a(i)*d(:,i);
g1=[g1 subs(subs(g,x1,X(1,i+1)),x2,X(2,i+1))];
m1=m2;m2=norm(g1(:,i+1));
k(i)=(m2/m1)^2;
i=i+1;
end
f1=subs(subs(f,x1,X(1,i)),x2,X(2,i));
elseif n==3
format long
syms x1 x2 x3
f=1/2*[x1,x2,x3]*G*[x1;x2;x3]+b'*[x1;x2;x3]+c;
g=[diff(f,x1);diff(f,x2);diff(f,x3)];
g1=subs(subs(subs(g,x1,X(1,1)),x2,X(2,1)),x3,X(3,1));
d=-g1;
a=-((X(:,1)'*G*d+b'*d)/(d'*G*d));
X(:,2)=X(:,1)+a*d;
g1=[g1 subs(subs(subs(g,x1,X(1,2)),x2,X(2,2)),x3,X(3,2))];
k=(norm(g1(:,2))/norm(g1(:,1)))^2;
m=norm(g1(:,2));
i=2;
while m>=e
d(:,i)=-g1(:,i)+k*d(:,i-1);
a(i)=-((X(:,i)'*G*d(:,i)+b'*d(:,i))/(d(:,i)'*G*d(:,i)));
X(:,i+1)=X(:,i)+a(i)*d(:,i);
g1=[g1 subs(subs(subs(g,x1,X(1,i+1)),x2,X(2,i+1)),x3,X(3,i+1))];
k=(norm(g1(:,i+1))/norm(g1(:,i)))^2;
m=norm(g1(:,i+1));
i=i+1;
end
f1=subs(subs(subs(f,x1,X(1,i)),x2,X(2,i)),x3,X(3,i));
end 参考技术A MATLAB中用遗传算法求解约束非线性规划问题 Solution of optimization with nonliear constraints programming by genetic alogorithm in MATLAB 作者:王勇, 期刊-核心期刊 哈尔滨商业大学学报（自然科学版）JOURNAL OF HARBIN UNIVERSITY OF COMMERCE(NATURAL SCIENCES EDITION) 2006年 第04期
- 约束优化问题的遗传算法求解 Genetic algorithm solution for constrained optimization 作者:宋松柏,蔡焕杰,康艳, 期刊-核心期刊 西北农林科技大学学报（自然科学版）JOURNAL OF NORTHWEST SCI-TECH UNIVERSITY OF AGRICULTURE AND FORESTRY(NATURAL本回答被提问者采纳 参考技术B 一般没有程序，你是不是问错了饿。啊