在有限元程序開發中，線性方程組的求解是一個重要組成部分。在百萬自由度大規模計算的情況下，線性方程組的高效快速求解對整個求解器的計算效率有著至關重要的作用。無論實際上計算的是線性問題，還是各種非線性問題，最終都需要落實到線性方程組的求解上去。非線性方程組的求解實際上往往就是多次求解線性方程組。

      目前，線性方程組的求解主要分為直接法和迭代法兩種。

 

      在之前的文章[數值算法與編程]高斯消去法中，我們討論的高斯消去法就是直接法的一種。而本文即將討論的共軛梯度法，是迭代法的一種，并且，其屬于目前求解對稱線性方程組的主要迭代方法。各大商業有限元軟件，在面臨對稱線性方程組的求解時幾乎都會選用各種變化形式的共軛梯度法進行求解。

      共軛梯度法的具體原理和算法如下：

      假定要求解的對稱線性方程組是：

     其中，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分解預處理等。

      以上，就是共軛梯度法求解線性方程組的一些簡要介紹，感謝您的閱讀！

【完】

 歡迎關注公眾號  有限元術