列主元高斯消元法¶

约 517 个字 85 行代码预计阅读时间 3 分钟

高斯消元过程中如果在消去过程中出现 0 主元或者是主元非常小的话，消去法将失败或者数值不稳定。这时可以采用选主元的方法，进行处理。

计算原理¶

设置增广矩阵 \(A_{p}=[A,b]\)
对 \(k=1\) 至 \(N-1\)，做3~6处理
设置一个元素 \(a_{max}=\vert A_{p}(k,k)\vert\)，及标号 \(ID_{max}=k\)
查找当列的最大元素，然后用标号 \(ID_{max}\) 记录下这个元素所在的行数
交换第 \(k\) 行与第 \(ID_{max}\) 行的所有数据。其他元素保持不变
完成5时已经完成了主元的选取，这时候可以按照上一节的消去法进行消去计算
完成以上步骤之后，已经形成上三角矩阵
对上三角矩阵进行回代，即可得到方程的解

graph TD
    A["01: 设置增广矩阵Ap=[A,b]"] --> B["02: k=1 to N-1"]
    B --> C["03: 初始化a_max=|Ap(k,k)|, IDmax=k"]
    C --> D["04: 查找第k列主元"]
    D --> E{"主元是否在k行?"}
    E -- 否 --> F["05: 交换行k与行IDmax"]
    E -- 是 --> G["06: 执行消去计算"]
    F --> G
    G --> H{"k < N-1?"}
    H -- 是 --> B
    H -- 否 --> I["07: 获得上三角矩阵"]
    I --> J["08: 回代求解"]

Fortran数值实现¶

subroutine m_gauss(a,b,x,n)
! 列主元高斯消元法求解线性方程组
    integer, intent(in) :: n
    real*8, intent(in) :: a(n,n),b(n)
    real*8, intent(out) :: x(n)

    integer :: i,k,j
    integer :: id_max ! 列主元对应的标号
    real*8 :: a_up(n,n),b_up(n)
    real*8 :: ab(n,n+1)
    real*8 :: vtemp1(n+1), vtemp2(n+1)
    ab(1:n,1:n) = a
    ab(:,n+1) = b

    !选主元
    do k = 1,n-1
        elmax = dabs(ab(k,k)) ! dabs()计算双精度浮点数的绝对值
        id_max = k
        ! 寻找最大元素对应的标号
        do i =k+1,n
            if (dabs(ab(i,k))>elmax) then
                elmax=Ab(i,k)
                id_max = i
            end if
        end do
        ! 与第k行进行交换
        vtemp1 = ab(k,:)
        vtemp2 = ab(id_max,:)

        ab(k,:) = vtemp2
        ab(id_max,:) = vtemp1
        ! 交换完毕，进行消元
        do i = k+1,n
            temp = ab(i,k)/ab(k,k)
            ab(i,:) = ab(i,:) - temp*ab(k,:)
        end do
    end do

    ! 打印ab，此时已经变换为上三角矩阵
    print *, 'Augmented matrix [A|b] after elimination:'
    do i = 1, n
        print '(7F12.6)', (ab(i,j), j = 1, n+1)
    end do

    a_up(:,:) = ab(:,1:n)
    b_up(:) = ab(:,n+1)

    call uptri(a_up,b_up,x,n)

end subroutine m_gauss

数值案例¶

\[ \mathbf { A } = \left( \begin{array} { r r r r r } { - 3 . 3 4 3 5 } & { - 0 . 4 9 4 6 } & { 0 . 3 8 3 4 } & { - 4 . 9 5 3 7 } & { - 2 . 4 0 1 3 } & { - 3 . 5 4 4 6 } \\ { 1 . 0 1 9 8 } & { - 4 . 1 6 1 8 } & { 4 . 9 6 1 3 } & { 2 . 7 4 9 1 } & { 3 . 0 0 0 7 } & { - 3 . 6 3 9 3 } \\ { - 2 . 3 7 0 3 } & { - 2 . 7 1 0 2 } & { 4 . 2 1 8 2 } & { 3 . 1 7 3 0 } & { - 0 . 6 8 5 9 } & { 3 . 6 9 2 9 } \\ { 1 . 5 4 0 8 } & { 4 . 1 3 3 4 } & { - 0 . 5 7 3 2 } & { 3 . 6 8 6 9 } & { 4 . 1 0 6 5 } & { 0 . 7 9 7 0 } \\ { 1 . 8 9 2 1 } & { - 3 . 4 7 6 2 } & { - 3 . 9 3 3 5 } & { - 4 . 1 5 5 6 } & { - 3 . 1 8 1 5 } & { 0 . 4 9 8 6 } \\ { 2 . 4 8 1 5 } & { 3 . 2 5 8 2 } & { 4 . 6 1 9 0 } & { - 1 . 0 0 2 2 } & { - 2 . 3 6 2 0 } & { - 3 . 5 5 0 5 } \end{array} \right) \]

\[ \mathbf { b } = \left( \begin{array} { c } { 9 . 0 5 3 7 } \\ { 0 . 0 2 2 8 } \\ { - 8 . 4 1 7 7 } \\ { - 4 . 6 3 8 0 } \\ { 1 0 . 5 5 7 5 } \\ { 9 . 8 2 5 2 } \end{array} \right) \, \]

主程序：

program main
    implicit none
    integer, parameter :: n = 6
    real*8 :: a(n,n), b(n), x(n)
    integer :: i

    ! 初始化矩阵A和向量b
    a = reshape([ &
        -3.3435d0, -0.4946d0,  0.3834d0, -4.9537d0, -2.4013d0, -3.5446d0, &
         1.0198d0, -4.1618d0,  4.9613d0,  2.7491d0,  3.0007d0, -3.6393d0, &
        -2.3703d0, -2.7102d0, -4.2182d0,  3.1730d0, -0.6859d0,  3.6929d0, &
         1.5408d0,  4.1334d0, -0.5732d0,  3.6869d0,  4.1065d0,  0.7970d0, &
         1.8921d0, -3.4762d0, -3.9335d0, -4.1556d0, -3.1815d0,  0.4986d0, &
         2.4815d0,  3.2582d0,  4.6190d0, -1.0022d0, -2.3620d0, -3.5505d0  &
    ], shape(a), order=[2,1])

    b = [9.0537d0,  0.0228d0, -8.4177d0, -4.6380d0, 10.5575d0, 9.8252d0]

    ! 调用高斯消元法求解
    call m_gauss(a, b, x, n)

    print *, "x:"
    do i = 1, n
        write(*, 100) "x(", i, ") = ", x(i)
    end do

100 format(1x, A, I1, A, F10.6)
end program main

输出：

 x:
 x(1) =   1.265159
 x(2) =   0.110264
 x(3) =  -1.245209
 x(4) =  -0.433799
 x(5) =  -0.990933
 x(6) =  -2.620118

注意：子程序中默认：implicit real*8(a-z)

参考文献¶

宋叶志,茅永兴,赵秀杰.Fortran 95/2003科学计算与工程[M].清华大学出版社,2011.