有限元方法中的数值技术：三角矩阵求解

从今天开始，木木将为大家分享自己在学习数值分析这门课过程中练习的Fortran代码，主要分享：

1. 线性方程组直接法求解
2. 线性方程组迭代法求解
3. 非线性方程求解
4. 非线性方程组求解
5. 数值积分
6. 矩阵特征值、特征向量
7. 插值、拟合

以上的内容均是有限元中常常用到的数值思想，可能大多数情况下就是调用一些前人写好的库，直接使用进行，不用关心内部实现原理。

但对于有些问题的思考、优化，如果脑袋里面没有一些数值技术的支撑，很难突破一些固有体系，学习更高阶的技能时也会非常受限。

最近，木木又重新捡起了Fortran语言，经典的公式语言，使用一段时间后才发觉是多么的优美，对于公式的每一处细节Fortran总能为你一一体现。

于是乎就用它进行后续一系列的编程，包括数值分析、结构动力学等研究生基础课程。

本系列推文非常适合研一、研二对数值计算感兴趣的同学，如果以后的研究领域与数值相关，非常建议通过编程来理解相关的核心课程，如果只是泛读文献、使用专业软件进行纯模拟仿真，深度是远远不够的，希望同学们可以在学习数值思想的同时也可以提升对语言的掌握程度。

本次分享的是三角方程组的解法，后续的很多解法都是在这个的基础上进行求解。

什么是三角方程组

三角方程组是指其系数矩阵呈三角形结构的方程组。根据非零元素分布位置不同，可分为：

1. 上三角方程组（Upper-Triangular System）
   系数矩阵        满足

      

      

   其一般形式为

      

2. 下三角方程组（Lower-Triangular System） 系数矩阵满足      ，其一般形式为

      

上三角形方程

计算原理

数值案例

subroutine uptri(a,b,x,n)
! 上三角方程组求解
    ! 输入参数：
    ! a(n,n): 上三角矩阵
    ! b(n): 右端向量
    ! n: 矩阵维度
    ! 输出参数：
    ! x(n): 解向量
    integer, intent(in) :: n
    real*8, intent(in) :: a(n,n),b(n)
    real*8, intent(out) :: x(n)

    integer :: i,k

    x(n) = b(n)/a(n,n)

    do i = n-1,1,-1
        x(i) = b(i)
        do k = i+1,n
            x(i) = x(i) - a(i,k)*x(k)
        enddo
        x(i) = x(i)/a(i,i)
    enddo
endsubroutine uptri

下三角方程

计算原理

数值案例

subroutine lowtri(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

    x(1) = b(1)/a(1,1)

    do k = 2,n
        x(k) = b(k)
        do i = 1,k-1
            x(k) = x(k) - a(k,i)*x(i)
        enddo
        x(k) = x(k)/a(k,k)
    enddo

endsubroutine lowtri

主程序调用：

program main
    use matrix
    implicitnone
    integer, parameter :: n = 4
    real*8 :: a_upper(n,n), a_lower(n,n), b_upper(n), b_lower(n), x(n), y(n)
    integer :: i
    
    ! 初始化上三角矩阵和右端向量
    a_upper = reshape([2.0, 0.0, 0.0, 0.0, &
                       1.0, 3.0, 0.0, 0.0, &
                       3.0, 2.0, 7.0, 0.0, &
                       4.0, 5.0, 3.0, 2.0], [n,n])
    b_upper = [50.0, 49.0, 53.0, 12.0]

    ! 初始化下三角矩阵和右端向量
    a_lower = reshape([2.0, 1.0, 1.0, 3.0, &
                       0.0, 3.0, 5.0, 2.0, &
                       0.0, 0.0, 1.0, 1.0, &
                       0.0, 0.0, 0.0, 4.0], [n,n])
    b_lower = [10.0, 23.0, 39.0, 43.0]


    ! 求解上三角方程组
    call uptri(a_upper, b_upper, x, n)
    ! 求解下三角方程组
    call lowtri(a_lower, b_lower, y, n)

    print *, "Solution of upper triangular system x:"
    do i = 1, n
        write(*, 100) "x(", i, ") = ", x(i)
    enddo
    
    print *, "----------------------------------------"
    
    print *, "Solution of lower triangular system y:"
    do i = 1, n
        write(*, 100) "y(", i, ") = ", y(i)
    enddo

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

输出：

 Solution of upper triangular system x:
 x(1) =   4.000000
 x(2) =   3.000000
 x(3) =   5.000000
 x(4) =   6.000000
 ----------------------------------------
 Solution of lower triangular system y:
 y(1) =   5.000000
 y(2) =   6.000000
 y(3) =   4.000000
 y(4) =   3.000000

 

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

参考文献

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