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

计算原理

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

$$\left(\begin{array}{llll}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& & \ddots & \\ a_{n1}& a_{n2}& \cdots & a_{nn}\end{array}\right)=\left(\begin{array}{llll}{l}_{11}& & & \\ l_{21}& l_{22}& & \\ \cdots & & \ddots & \\ l_{n1}& l_{n2}& \cdots & l_{nn}\end{array}\right)\left(\begin{array}{llll}{l}_{11}& l_{21}& \cdots & l_{n1}\\ & l_{22}& \cdots & l_{n2}\\ & & \ddots & \\ & & & l_{nn}\end{array}\right)$$

1. $l_{11}=\sqrt{a_{11}}$
2. $l_{i1}=\frac{a_{i1}}{l_{11}},i=2,N$
3. 对 $j=2,N$，做A~B处理：

A:
$$l_{jj}=\sqrt{a_{jj}-\sum _{k=1}^{j-1}{l}_{jk}^2}$$

B:
$$l_{ij}=\frac{\left(a_{ij}-{\sum }_{k=1}^{j-1}{l}_{ik}{l}_{jk}\right)}{l_{jj}},i=j+1,N$$

Fortran数值实现

subroutine cholesky(a,l,n)
! cholesky矩阵分解(只适用于对称正定矩阵)! A = L*L^Tinteger, intent(in) :: nreal*8, intent(in) :: a(n,n)real*8, intent(out) :: l(n,n)integer :: i,j,kl = 0l(1,1) = sqrt(a(1,1))l(2:,1) = a(2:,1)/l(1,1)do j = 2,ns = 0do k = 1,j-1s = s + l(j,k)**2end dol(j,j) = sqrt(a(j,j)-s)do i = j+1,ns = 0do k = 1,j-1s = s + l(i,k) * l(j,k)end dol(i,j) = (a(i,j)-s)/l(j,j)end doend do
end subroutine cholesky

数值案例

$$\mathbf{A}=\left(\begin{array}{cccc}233.4615& 113.8423& 256.0623& 145.0697\\ 113.8423& 78.6033& 127.4298& 95.3089\\ 256.0623& 127.4298& 281.4721& 164.8676\\ 145.0697& 95.3089& 164.8676& 181.2339\end{array}\right)$$

主程序：

program mainimplicit noneinteger,  parameter :: n = 4real*8 :: a(n,n),l(n,n)integer :: i, ja = 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, nwrite(*, 100) (l(i,j), j = 1, n)end do100 format(4F12.6)
end program main

输出：

 Lower Triangular Matrix L from Cholesky decomposition:15.279447    0.000000    0.000000    0.0000007.450682    4.805272    0.000000    0.00000016.758610    0.534147    0.579450    0.0000009.494434    5.112904    5.217081    6.142468

求解矩阵方程

对于$\mathbf{Ax}=\mathbf{b}$，进行cholesky分解：

$$\mathbf{L}{\mathbf{L}}^T\mathbf{x}=\mathbf{b}，等价于\{\begin{array}{l}\mathbf{Ly}=\mathbf{b}\\ {\mathbf{L}}^{\mathrm{T}}\mathbf{x}=\mathbf{y}\end{array}$$

对于：

$$\mathbf{A}=\left(\begin{array}{lll}4& -1& 1\\ -1& 4.25& 2.75\\ 1& 2.75& 3.5\end{array}\right),\mathbf{b}=\left(\begin{array}{l}4\\ 6\\ 7.25\end{array}\right)$$

program mainimplicit noneinteger,  parameter :: n = 3real*8 :: a(n,n), l(n,n), b(n), x(n), y(n)integer :: ia = 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 = bcall lowtri(l, b, y, n)! 解 L^T x = ycall uptri(transpose(l), y, x, n)print *, "Solution vector x:"write(*, 100) x100 format(3F12.6)
end program main

输出：

 Solution vector x:1.000000    1.000000    1.000000

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

参考文献

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