Doolittle矩阵分解¶
约 450 个字 102 行代码 预计阅读时间 3 分钟
算法原理¶
L、U分解后：
\[
\mathbf{L}=\left( \begin{array}{l}
1& 0& 0& 0\\
l_{21}& 1& 0& 0\\
\vdots& ...& \ddots& 0\\
l_{n1}& l_{n2}& \cdots& 1\\
\end{array} \right)
\]
\[
\mathbf{U}=\left( \begin{array}{l}
u_{11}& u_{12}& \cdots& u_{1n}\\
0& u_{22}& \cdots& u_{2n}\\
0& 0& \ddots& \vdots\\
0& 0& 0& u_{nn}\\
\end{array} \right)
\]
- 计算\(LU\)分解中的\(U\)的第一行元素和\(L\)的第一列
\[
\begin{array} { l } { { u _ { 1 i } = a _ { 1 i } , i = 1 , \cdots , n } } \\ { { l _ { j 1 } = \displaystyle \frac { a _ { j 1 } } { u _ { 1 1 } } , j = 1 , \cdots , n } } \end{array}
\]
- 对于\(k=2,···,n\)计算以及\(U\)的第\(k\)行元素：
\[
u _ { k j } = a _ { k j } - \sum _ { m = 1 } ^ { k - 1 } l _ { k m } u _ { m j } \, , \quad j = k , k + 1 , \cdots , n
\]
- 以及\(L\)的第\(k\)列元素：
\[
l _ { i k } = \frac { a _ { i k } - \sum _ { m = 1 } ^ { k - 1 } l _ { i m } u _ { m k } } { u _ { k k } } , \quad i = k + 1 , \cdots , n
\]
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)
end do
u(k,j) = a(k,j) - s
end do
do i = k+1,n
s = 0
do m = 1,k-1
s = s + l(i,m) * u(m,k)
end do
l(i,k) = (a(i,k)-s)/u(k,k)
end do
end do
end subroutine doolittle
数值案例¶
\[
\begin{array} { r } { \mathbf { A } = \left( \begin{array} { c c c } { 6 } & { 3 } & { - 8 } \\ { 1 5 } & { 5 } & { 2 } \\ { 2 } & { 2 } & { 7 } \end{array} \right) } \end{array}
\]
主程序：
program main
implicit none
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)
end do
print *
call doolittle(a, l, u, n)
print *, "Lower Triangular Matrix L:"
do i = 1, n
write(*, 100) (l(i,j), j = 1, n)
end do
print *
print *, "Upper Triangular Matrix U:"
do i = 1, n
write(*, 100) (u(i,j), j = 1, n)
end do
100 format(3F12.6)
end program 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
求解矩阵方程¶
对于对于\(\mathbf { A x } = \mathbf { b }\)，进行Doolittle分解：\(\mathbf { L U x } = \mathbf { b }\)
分两步求解：
- \(\mathbf { L } \mathbf { y } = \mathbf { b }\)
- \(\mathbf { U } \mathbf { x } = \mathbf { y }\)
\[
\mathbf { A } = \left( \begin{array} { c c c c } { 6 . 5 5 7 4 } & { 6 . 7 8 7 4 } & { 6 . 5 5 4 8 } & { 2 . 7 6 9 2 } \\ { 0 . 3 5 7 1 } & { 7 . 5 7 7 4 } & { 1 . 7 1 1 9 } & { 0 . 4 6 1 7 } \\ { 8 . 4 9 1 3 } & { 7 . 4 3 1 3 } & { 7 . 0 6 0 5 } & { 0 . 9 7 1 3 } \\ { 9 . 3 3 9 9 } & { 3 . 9 2 2 3 } & { 0 . 3 1 8 3 } & { 8 . 2 3 4 6 } \end{array} \right) , \mathbf { b } = \left( \begin{array} { c } { 1 3 0 . 3 2 4 2 } \\ { 4 2 . 9 3 4 8 } \\ { 1 4 9 . 9 8 9 3 } \\ { 8 3 . 1 9 5 3 } \end{array} \right)
\]
program main
implicit none
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
100 format(4F12.6)
end program main
输出：
注意：子程序中默认：
implicit real*8(a-z)
参考文献¶
- 宋叶志,茅永兴,赵秀杰.Fortran 95/2003科学计算与工程[M].清华大学出版社,2011.