有限元方法中的数值技术：高斯消元法

计算原理

1. 将系数矩阵与右端向量组合为增广矩阵      
2. 进行初等行/列变换，直至变换为上三角矩阵
3. 将上一节的上三角矩阵求解方法代入即可

第    次初等行变换中，令    ，第    行的消去公式为：

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

数值案例

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,:)
        enddo
    enddo
    a_up(:,:) = ab(:,1:n)
    b_up(:) = ab(:,n+1)

    call uptri(a_up,b_up,x,n)
endsubroutine gauss

主程序：

program main
    use matrix
    implicitnone
    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)
    enddo

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

输出：

 x:
 x(1) =  -3.378184
 x(2) =   2.942842
 x(3) =  -1.887843
 x(4) =   0.285336

 

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

参考文献

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