【數值算法】系數矩陣非對稱時，線性方程組如何求解？-穩定雙共軛梯度法(Bicgstab)求解線性方程組

寒江雪_123

2022年6月30日 20:32

瀏覽：2982 評論：12 收藏：2

在前面的文章和中表明共軛梯度法是求解對稱正定線性方程組的一種有效方法，當針對不同的系數矩陣采用不同的預處理方式時，其可以以較少的迭代次數獲得較高精度的解。然而，該方法的一個缺點就是其只能適用于對稱正定系數矩陣，當系數矩陣不再是對稱正定時，此方法可能失效。

以下舉例：

【數值算法】系數矩陣非對稱時，線性方程組如何求解？-穩定雙共軛梯度法(Bicgstab)求解線性方程組的圖2

上面矩陣A為非對稱矩陣，采用共軛梯度法求解過程如下：

【數值算法】系數矩陣非對稱時，線性方程組如何求解？-穩定雙共軛梯度法(Bicgstab)求解線性方程組的圖5

該方程組采用共軛梯度法迭代4862次依然未收斂。因此，對于該非對稱方程，可以認為，共軛梯度法幾乎是失效的。

在實際工程中，有限元方法形成的剛度系數以對稱正定居多，但是實際上也存在非對稱的可能，例如，當材料本構采用摩爾-庫倫本構時，其形成的剛度矩陣就有可能會是非對稱的，此時如果是使用商業軟件，應當在軟件中選擇非對稱求解器。如果是自主編程且采用迭代法求解線性方程組，則需要找到一種適用于非對稱矩陣的求解方法。

【數值算法】系數矩陣非對稱時，線性方程組如何求解？-穩定雙共軛梯度法(Bicgstab)求解線性方程組的圖7

常見的非對稱系數矩陣求解方法主要有：廣義最小殘差法(GMRES)，雙共軛梯度法（Bicg）穩定雙共軛梯度法(BiCGStab)，穩定混合雙共軛梯度法（BiCGStab(l)），這些方法相對于常規的共軛梯度法在推導上均增加了一些難度，實際推導往往較為復雜。本文不展開推導，僅對穩定雙共軛梯度法(BiCGStab)的偽代碼作簡要粘貼。

BiCGStab法的具體計算過程如下：

【數值算法】系數矩陣非對稱時，線性方程組如何求解？-穩定雙共軛梯度法(Bicgstab)求解線性方程組的圖9

具體代碼：

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

依據上述過程編寫程序，計算前述非對稱矩陣線性方程組求解結果：