高斯消元法¶

约 454 个字 61 行代码预计阅读时间 2 分钟

计算原理¶

将系数矩阵与右端向量组合为增广矩阵 \([A|b]\)
进行初等行/列变换，直至变换为上三角矩阵
将上一节的上三角矩阵求解方法代入即可

\[ \left[ \mathbf { A } ^ { ( n ) } \quad \mathbf { b } ^ { ( n ) } \right] = \left( \begin{array} { c c c c c } { a _ { 1 1 } ^ { ( 1 ) } } & { a _ { 1 2 } ^ { ( 1 ) } } & { \cdots } & { a _ { 1 n } ^ { ( 1 ) } } & { b _ { 1 } ^ { ( 1 ) } } \\ { 0 } & { a _ { 2 2 } ^ { ( 2 ) } } & { \cdots } & { a _ { 2 n } ^ { ( 2 ) } } & { b _ { 2 } ^ { ( 2 ) } } \\ { \vdots } & { \vdots } & { \ddots } & { \vdots } & { \vdots } \\ { 0 } & { 0 } & { \cdots } & { a _ { n n } ^ { ( n ) } } & { b _ { n } ^ { ( n ) } } \end{array} \right) \]

第\(k\)次初等行变换中，令\(l_{ik}=\frac{a_{ik}^k}{a_{kk}^k}\)，第\(i\)行的消去公式为：

\[ \left\{ \begin{array} { l l } { a _ { i j } ^ { ( k + 1 ) } = a _ { i j } ^ { ( k ) } - l _ { i k } a _ { i j } ^ { ( k ) } \quad i , j = k + 1 , k + 2 , . . . , n } \\ { b _ { i } ^ { k + 1 } = b _ { i } ^ { ( k ) } - l _ { i k } b _ { i } ^ { ( k ) } \quad i = k + 1 , k + 2 , . . . , n } \end{array} \right. \]

当转换为上三角矩阵后，再使用回代公式：

\[ \left\{ \begin{array} { l l } { \displaystyle x _ { n } = \frac { b _ { n } } { a _ { n n } } } \\ { \displaystyle x _ { i } = \frac { b _ { i } - \sum _ { k = i + 1 } ^ { n } a _ { i k } x _ { k } } { a _ { i i } } , i = n - 1 , n - 2 , \cdots , 1 } \end{array} \right. \]

数值案例¶

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

subroutine gauss(a,b,x,n)
! 高斯消元法求解线性方程组
! 输入参数：
! a(n,n): 系数矩阵
! b(n): 右端向量
! n: 矩阵维度
! 输出参数：
! x(n): 解向量
    integer :: i,k,n
    real*8 :: a(n,n),b(n),x(n)
    real*8 :: a_up(n,n),b_up(n)
    ! 增广矩阵 [A|b]
    real*8 :: ab(n,n+1)

    ab(1:n,1:n) = a
    ab(:,n+1) = b

    ! 增广矩阵行变换
    do k = 1,n-1
        do i = k+1,n
            temp = ab(i,k)/ab(k,k)
            ab(i,:) = ab(i,:) - temp*ab(k,:)
        end do
    end do
    a_up(:,:) = ab(:,1:n)
    b_up(:) = ab(:,n+1)

    call uptri(a_up,b_up,x,n)
end subroutine gauss

主程序：

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

    ! 初始化矩阵A和向量b
    a = reshape([ &
        -13.9211d0, 18.3862d0, -4.1683d0, -6.0438d0, &
        21.7540d0, -26.0893d0, 3.9325d0, 6.7018d0, &
        -14.8743d0, -5.6866d0, -33.3169d0, -32.9591d0, &
        -7.9025d0, 4.4451d0, 41.7098d0, -23.3378d0], &
        [n, n])

    b = [136.8721d0, -126.8849d0, 100.4524d0, 95.7019d0]

    ! 调用高斯消元法求解
    call 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) =  -3.378184
 x(2) =   2.942842
 x(3) =  -1.887843
 x(4) =   0.285336

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

参考文献¶

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