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

计算原理

对于对称正定矩阵    ，存在于一个实的下三角矩阵    ，使得    ，若限定    的对角元素为正时，这种分解是唯一的。

1.      
2.      
3. 对       ，做A~B处理：

A:

B:

Fortran数值实现

subroutine cholesky(a,l,n)
! cholesky矩阵分解(只适用于对称正定矩阵)
! A = L*L^T
integer, intent(in) :: n
real*8, intent(in) :: a(n,n)
real*8, intent(out) :: l(n,n)

integer :: i,j,k

 l = 0
 l(1,1) = sqrt(a(1,1))
 l(2:,1) = a(2:,1)/l(1,1)

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

数值案例

主程序：

program main
    implicitnone
    integer,  parameter :: n = 4
    real*8 :: a(n,n),l(n,n)
    integer :: i, j
    
    a = reshape([ &
        233.4615d0, 113.8423d0, 256.0623d0, 145.0697d0, &
        113.8423d0,  78.6033d0, 127.4298d0,  95.3089d0, &
        256.0623d0, 127.4298d0, 281.4721d0, 164.8676d0, &
        145.0697d0,  95.3089d0, 164.8676d0, 181.2339d0  &
    ], shape(a), order=[2,1])
    
    call cholesky(a, l, n)
    
    print *, "Lower Triangular Matrix L from Cholesky decomposition:"
    do i = 1, n
        write(*, 100) (l(i,j), j = 1, n)
    enddo
    
100format(4F12.6)
endprogram main

输出：

 Lower Triangular Matrix L from Cholesky decomposition:
   15.279447    0.000000    0.000000    0.000000
    7.450682    4.805272    0.000000    0.000000
   16.758610    0.534147    0.579450    0.000000
    9.494434    5.112904    5.217081    6.142468

求解矩阵方程

对于    ，进行cholesky分解：

对于：

program main
    implicitnone
    integer,  parameter :: n = 3
    real*8 :: a(n,n), l(n,n), b(n), x(n), y(n)
    integer :: i
    
    a = reshape([ &
        4.0d0,   -1.0d0,   1.0d0,  &
        -1.0d0,  4.25d0,  2.75d0, &
        1.0d0,   2.75d0,  3.5d0   &
    ], shape(a), order=[2,1])
    
    b = [4.0d0, 6.0d0, 7.25d0]

    call cholesky(a, l, n)
    ! 解 Ly = b
    call lowtri(l, b, y, n)
    ! 解 L^T x = y
    call uptri(transpose(l), y, x, n)
    
    print *, "Solution vector x:"
    write(*, 100) x
    
100format(3F12.6)
endprogram main

输出：

 Solution vector x:
    1.000000    1.000000    1.000000