有限元方法中的数值技术：Doolittle矩阵分解

算法原理

L、U分解后：

• 计算      分解中的      的第一行元素和      的第一列

• 对于      计算以及      的第      行元素：

• 以及      的第      列元素：

Fortran数值实现

subroutine doolittle(a,l,u,n)
! doolittle矩阵分解
! A = LU
integer, intent(in) :: n
real*8, intent(in) :: a(n,n)
real*8, intent(out) :: l(n,n),u(n,n)

integer :: i,j,k,m

! U的第一行和L的第一列
 u(1,:) = a(1,:)
 l(:,1) = a(:,1)/u(1,1)

do k = 2,n
  l(k,k) = 1
do j = k,n
   s = 0
   do m = 1,k-1
    s = s + l(k,m) * u(m,j)
   enddo
   u(k,j) = a(k,j) - s
enddo
do i = k+1,n
   s = 0
   do m = 1,k-1
    s = s + l(i,m) * u(m,k)
   enddo
   l(i,k) = (a(i,k)-s)/u(k,k)
enddo
enddo
endsubroutine doolittle

数值案例

主程序：

program main
    implicitnone
    integer, parameter :: n = 3
    real*8 :: a(n,n), l(n,n), u(n,n)
    integer :: i, j
    
    a = reshape([ 6.0d0, 3.0d0, -8.0d0, &
                  15.0d0, 5.0d0, 2.0d0, &
                  2.0d0, 0.0d0, 7.0d0  &
                ], shape(a), order=[2,1])
    
    print *, "Original Matrix A:"
    do i = 1, n
        write(*, 100) (a(i,j), j = 1, n)
    enddo
    print *
    
    call doolittle(a, l, u, n)
    
    print *, "Lower Triangular Matrix L:"
    do i = 1, n
        write(*, 100) (l(i,j), j = 1, n)
    enddo
    print *
    
    print *, "Upper Triangular Matrix U:"
    do i = 1, n
        write(*, 100) (u(i,j), j = 1, n)
    enddo
    
100format(3F12.6)
endprogram main

输出：

 Lower Triangular Matrix L:
    1.000000    0.000000    0.000000
    2.500000    1.000000    0.000000
    0.333333    0.400000    1.000000

 Upper Triangular Matrix U:
    6.000000    3.000000   -8.000000
    0.000000   -2.500000   22.000000
    0.000000    0.000000    0.866667

求解矩阵方程

对于对于    ，进行Doolittle分解：    

分两步求解：

program main
    use matrix
    implicitnone
    integer,  parameter :: n = 4
    real*8 :: a(n,n), l(n,n), u(n,n), b(n), x(n), y(n)
    integer :: i
    
    a = reshape([ &
        6.5574d0, 6.7874d0, 6.5548d0, 2.7692d0, &
        0.3571d0, 7.5774d0, 1.7119d0, 0.4617d0, &
        8.4913d0, 7.4313d0, 7.0605d0, 0.9713d0, &
        9.3399d0, 3.9223d0, 0.3183d0, 8.2346d0  &
    ], shape(a), order=[2,1])
    
    b = [130.3242d0, 42.9348d0, 149.9893d0, 83.1953d0]

    call doolittle(a, l, u, n)

    ! 解 Ly = b
    call lowtri(l, b, y, n)

    ! 解 U x = y
    call uptri(u, y, x, n)
    
    print *, "Solution vector x:"
    write(*, 100) x
    
100format(4F12.6)
endprogram main

输出：

 Solution vector x:
    6.948332    3.170983    9.502135    0.344460

 

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

参考文献

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