在之前的文章[數值算法與編程]高斯消去法 中，本公眾號編寫了高斯消去法求解線性方程組的具體代碼。其具體算法如下：

（1）消元部分

（2）回代部分

 

   很明顯，對于某些矩陣，使用上述算法可能會出現a(i,i)為0的情況，而一旦出現這種情況，該算法實際上就無法繼續(xù)進行求解。

   以以下方程組為例：

上述方程組的系數矩陣為：

第一次消元后，系數矩陣變?yōu)?/span>

   顯然，由于A1(2,2)=0，下一次消元已經無法進行，因此直接采用高斯消去法對該方程組是無法進行求解的。

 

   針對該問題，改進的算法叫選主元高斯消去法。

具體步驟如下：

1.設置增廣矩陣AB=[A,B]

2.對增廣矩陣的第i（i=1~N-1）列進行以下處理：

2.1 設amax=ABS(AB(i,i))，IDmax=i

2.2對第i列的對角線以下的元素進行遍歷，如果有元素的絕對值大于對角線元素的絕對值（amax），則獲得該元素所在行并將其行數賦值給IDMAX，并將AB矩陣中該行（即IDMAX行）與第K行進行元素交換。

對變換后的矩陣按照普通的高斯消去法進行求解。

具體代碼：

 program Console1_xuanzhuyuangaosi

 

    implicit none

    integer,parameter::n=5

    real::a(n,n)

    real::b(n)

    integer::i

    real::x(n)

 

 

 

 

    open(10,file='A.txt')

    read(10,*)(a(i,:),i=1,n)

 

    open(12,file='B.txt')

    read(12,*)(b(i),i=1,n)

 

    write(*,*)(a(i,:),i=1,n)

    write(*,*)(b(i),i=1,n)

    ! Variables

 

    call select_gauss(a,b,x,n)

end program Console1_xuanzhuyuangaosi

 

其中A.txt中的內容：

1     2     3     4     5

2     4     5     3     6

1     1     1     1     1

2     3     6     4     2

3     3     2     4     2

 

  B.txt中的內容：

8

8

8

8

8

  

 

 

 

 

    subroutine select_gauss(A,B,X,N)

    implicit none

    real::A(n,n)

    real::B(n)

    real::X(n)

    real::Ap(n,n+1)

    integer::k,IDmax,i,j

    real::amax

    real::temp1(n+1),temp2(n+1)

    real::colmax

    real::xishu,sum

    integer::n

    Ap(:,1:n)=A(:,:)

    Ap(:,n+1)=B(:)

   

 

  

    write(*,"(7f20.6)")(ap(i,:),i=1,n)

 

    do k=1,n-1

        amax=abs(Ap(k,k))

        IDmax=k

        do i=k+1,n

            if(abs(Ap(i,k))>amax) then

                IDmax=i

                amax=ap(i,k)

            endif

        enddo 

        temp1(:)=Ap(k,:)

        temp2(:)=ap(idmax,:)

        Ap(k,:)=temp2(:)

        Ap(IDmax,:)=temp1(:) 

        do i=k+1,n

            xishu=ap(i,k)/ap(k,k)

            ap(i,:)=ap(i,:)-ap(k,:)*xishu

        end do

    end do

    

     write(*,*)"交換后"

        write(*,"(7f12.4)")(ap(i,:),i=1,n)

 

    x(n)=ap(n,n+1)/ap(n,n)

    

    do i=n-1,1,-1    

        sum=0.0

        do j=i+1,n

            sum=sum+ap(i,j)*x(j)

        end do

        x(i)=(ap(i,n+1)-sum)/(ap(i,i))

    end do

 

    write(*,*)"最終的解是：",(x(i),i=1,n)

end subroutine select_gauss

 

求解結果：

【完】