【数值算法】共轭梯度法求解线性方程组

   

      在有限元程序开发中，线性方程组的求解是一个重要组成部分。在百万自由度大规模计算的情况下，线性方程组的高效快速求解对整个求解器的计算效率有着至关重要的作用。无论实际上计算的是线性问题，还是各种非线性问题，最终都需要落实到线性方程组的求解上去。非线性方程组的求解实际上往往就是多次求解线性方程组。

      目前，线性方程组的求解主要分为直接法和迭代法两种。

 

      在之前的文章[数值算法与编程]高斯消去法中，我们讨论的高斯消去法就是直接法的一种。而本文即将讨论的共轭梯度法，是迭代法的一种，并且，其属于目前求解对称线性方程组的主要迭代方法。各大商业有限元软件，在面临对称线性方程组的求解时几乎都会选用各种变化形式的共轭梯度法进行求解。

      共轭梯度法的具体原理和算法如下：

      假定要求解的对称线性方程组是：

     其中，A是对称正定的系数矩阵。

     则实际上待求的解也是方程

      取得最小值的时候的解。

     求该方程的最小值的常见方法是最速下降法，该方法算法伪代码如下：

    

      该方法实际上是沿着负梯度方向进行搜索，直至残量接近0，较为简便，但是在条件数很大时，该方法收敛很慢。

      共轭梯度法，就是在上述的最速下降法的基础上进行改进的，其不沿着负梯度方向搜索，而沿着与梯度共轭的方向进行搜索，直至搜索置残量接近0，具体算法的伪代码如下：

 

具体代码：

program pcg_main

    implicit none

    integer,parameter::n=8

    real::a(n,n),b(n)

    real::x(n),x0(n)

    integer::i,j

    

    a=reshape((/6,0,1,2,0,0,2,1, &

                0,5,1,1,0,0,3,0,&

                1,1,6,1,2,0,1,2, &

                2,1,1,7,1,2,1,1,&

                0,0,2,1,6,0,2,1,&

                0,0,0,2,0,4,1,0,&

                2,3,1,1,2,1,5,1,&

                1,0,2,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 pcg(a,b,x,x0,n)

    end program pcg_main

subroutine pcg(a,b,x,x0,n)

    implicit none

    integer,intent(in)::n

    real,intent(in)::a(n,n),b(n),x0(n)

    real,intent(out)::x(n)

    real::xkp1(n)

    real::rk(n),r0(n),rkp1(n)

    real::ak(1),bk  

    real::pk(n),p0(n),pkp1(n)

    integer::k

    integer::count

    real::xk(n)

    real::dx=0.0

    integer::i

    real::tol=1.0e-8

    real::pk1(1,n)

   

    r0=b-matmul(a,x0)

 !   write(*,"(*(g0,1x))")"r0=",r0

    p0=r0

    pk=p0

    xk=x0

    rk=r0

    count=0

    ak=0.0

    xkp1=0.0

    pk1=0.0

 

    outer:do

        

        pk1(1,:)=pk(:)

        ak(:)=(norm2(rk)**2.0)/(matmul(pk1,matmul(a,pk)))    

        xkp1=xk+ak(1)*pk  

        write(*,*)count,xkp1

        dx=norm2(rk)**2.0    

        write(*,*)"dx=",dx

        if(dx<tol) then

            write(*,*)"converged,exit!"

            exit outer

        endif

 

        rkp1=rk-ak(1)*matmul(a,pk)

        bk=norm2(rkp1)**2.0/norm2(rk)**2.0

        pkp1=rkp1+bk*pk

        rk=rkp1

        xk=xkp1

        pk=pkp1

        count=count+1

        write(*,*)"第",count,"次迭代"

    enddo outer

    x=xkp1

    write(*,*)"最终解是",x

    end subroutine pcg

计算结果：

Matlab结果：

 

      可以看到，在容差设置为1.0e-8的情况下，经过9次迭代，获得了收敛解。实际上从数学上可以证明，如果对称正定矩阵A的阶数是n，通常最多第n+1次迭代就可获得收敛。

      共轭梯度法对于求解对称正定线性方程组十分有效，基本上是对称正定线性方程组迭代法求解的不二选择。在实际的大型线性方程组求解中，会进一步对该算法进行预处理改进，使得其迭代次数更少，求解更快，常见的预处理方法有对角线预处理，不完全Cholesky分解预处理等。

      以上，就是共轭梯度法求解线性方程组的一些简要介绍，感谢您的阅读！

【完】

 欢迎关注公众号  有限元术