[数值算法与编程]选主元高斯消去法

   在之前的文章[数值算法与编程]高斯消去法 中，本公众号编写了高斯消去法求解线性方程组的具体代码。其具体算法如下：

（1）消元部分

（2）回代部分

 

   很明显，对于某些矩阵，使用上述算法可能会出现a(i,i)为0的情况，而一旦出现这种情况，该算法实际上就无法继续进行求解。

   以以下方程组为例：

上述方程组的系数矩阵为：

第一次消元后，系数矩阵变为

   显然，由于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

 

求解结果：

【完】