数值分析第4讲矩阵的三角分解法

a2n ann
l21 ln1
1 ln2
u11 u12
u22
1
u1n u2n unn
比较第1行： a1 j u1 j 比较第1列：ai1 li1u11
j 1, , n i 2, , n
u1 j a1 j
li1
ai1 u11
数值分析
比较第2行： a2 j l21u1 j u2 j j 2, , n 比较第2列： ai2 li1u12 li2u 22 i 3, , n
k 1
(2) 选主元：取ir使得
sir
max
rin
si
，Ip(r)
ir
(3) 交换A的第r行与第ir行：ari air ,r , (i 1,2, , n)
(4) 计算U的第r行, L的第r列元素
arr urr sr
air lir si / urr air / arr，(i r 1, , n, r n)
r 1 及 urr arr lrkukr .
k 1
希望urr为大的数！
数值分析
为了避免用小的数urr作除数，引进量
r 1 si air likukr (i r, , n).
k 1
当不选主元时，
urr sr , lir si / sr (i r 1, , n).
当选主元时，取ir使得
r 1
这时lir 1
ari uri ari lrkuki (i r 1, , n, r n) k 1
数值分析
求解Ly=Pb及Ux=y的算法：
2. 对i 1,2, ,n 1,
(1) t Ip(i) (2) 如果i t则转(3)
bi bt
(3) 继续循环
i 1 3. bi bi lik bk (i 2,3, , n)
No pivoting
Directly factorizatio
n
column pivoting
Eliminating elements in a
row except for the diagonal element
A symmetric
and positive definite
No pivoting
d1 D1 0, di Di / Di1 0,(i 2, , n).
数值分析
d1
D
d1
dn
d1
dn
11
D2D2
dn
A
LDLT
LD
1 2
D
1 2
LT
1
(LD2
1
)(LD2
)T
L1LT1
定理11(对称正定矩阵的三角分解或Cholesky分解) 设A为 n阶对称正定矩阵, 则存在实的非奇异下三角矩阵L使得A
bi lij y j
yi
j 1
lii
, i 1, , n
n
xi yi uij x j , i n, ,1 j i 1
数值分析
数值分析
定理2.4 若矩阵A非奇异，则存在置换矩阵Q， 使得QA可以作Doolittle分解
QA=LU 其中L是单位下三角矩阵，U是上三角矩阵。
数值分析
选主元直接三角分解法
sir
max si ,交换A的第r行与
rin
第ir行,将sir 调到(r, r)位置，(i, j)处新元素仍记为lij和aij,
再进行第r步分解计算.
数值分析
选主元三角分解算法：PA A, 整型Ip(n)记录主行, x b.
1. 对r 1,2, ,n,
(1) 计算si :
r 1 air si air likukr (i r, , n).
a1 b1
1
1 1
c2
a2
cn
bn1 an
2
2
n
n
1
n 1
1
ii
ci ai
, i 2,
ci i1 ,
,n i 1,
,n
i bi i , i 1, , n
, c1 0
yi
(
fi
ci yi1)
i
xi yi i xi1 ,
, i 1, , n
i n, ,1 (n 0)
j k, ,n
k 1
ukj akj lkrurj r 1
比较第k列：
k 1
aik lirurk likukk r 1
K-1次
k 1
aik lirurk
i k 1, , n lik
Leabharlann Baidur 1
ukk
K-1＋1次
数值分析
分解过程完毕，加上两次反代过程
i 1
yi bi lij y j j 1
数值分析
所以，有计算过程如下：
i ai ci i1
i
yi
bi
i
( fi
ci
yi 1 )
i
xk yk k xk1 ,
1 2 n和y1 y2 yn 为追；xn xn1 x1为赶.
i 1, , n
k n, ,1
说明: 稳定性（对角占优）；运算量5n-4次乘除法；
存贮（一维数组）.
GPBi-CG Bi-CGSTAB(L)
Bi-CGSTAB2 Bi-CGSTAB
CGS
GMRES
Conjugate Gradient method for symmetric systems
Bi-CG
数值分析
回顾高斯消去法
Gauss消元法的第k步：
第i行
a (k ) ik
第k行, i
k
1,
,n
b y
数值分析
k 1
比较第k列： aik lirurk lik r 1
i k, ,n
k 1
lik aik lirurk r 1
比较第k行：
k 1
akj lkrurj lkkukj r 1
j k 1, , n
k 1
akj lkrurj
ukj
r 1
lkk
两次反代过程
i 1
数值分析
平方根法
应用有限元法解结构力学问题时，最后归结为求解线性 代数方程组，系数矩阵往往对称正定。平方根法是一种 对称正定矩阵的三角分解法，广泛用于求解系数矩阵为 对称正定的线性代数方程组。
设A为对称矩阵,且顺序主子式不为零，则
u11
U
u22
1
u12 u11
1
unn
u23 u22
yi
bi
y1 b1
i 1
lij y j , i 2,L
j 1
,n
xn yn unn
xi
yi
n
uij y j
u ii , i n 1,L ,1
ji1
数值分析
1、Doolittle分解
L为单位下三角，U为上三角
a11 a12
a21 a22
an1
an2
a1n 1
ln2
u1n u2n unn
数值分析
例：用直角三角分解法解
1 2 3 x1 14
2
5
2
x2
18
3 1 5 x3 20
解：用分解计算公式得
1 0 0 1 2 3
A 2
1
0 0 1
4
LU .
3 5 1 0 0 24
求解 Ly (14,18, 20)T y (14, 10, 72)T ,
LLT ,当限定L对角元素为正时，这种分解是唯一的.
数值分析
l11 0 A l21 l22
ln1 ln2
0 l11 l21
0
l22
0
lnn 0 0
ln1
ln
2
,
lnn
Cholesky分解 : 对于j 1,2, ,n,
j 1 1
1.
l jj (a jj
l
2 jk
, i 1, , n
总运算量为：
n
yi uij x j
xi
j i1
uii
, i n, ,1
n1 (n k 1)(k 1) (n k)k n(n 1) n(n 1) n3 n2 n
k 1
2
2
3
3
数值分析
存储在矩阵的原来位置，且不影响计算
u11 u12
l21 u22
ln1
k 1 n
4. bn bn / unn , bi (bi uikbk ) / uii k i 1
(i n 1, ,1).
A1 U 1L1P
数值分析
数值分析
数值分析
三对角阵的追赶法（A的前n-1个顺序主子式非零）
在数值求解常微分方程边值问题、热传导方程和建 立三次样条函数时，都会要解三对角方程组：AX=b
3：对上半带宽为s、下半带宽为r的带状矩阵作LU分解，那么L 为下半带宽为r的下三角矩阵，U为上半带宽为s的上三角矩阵；
数值分析
直接法内容总结
4：对上半带宽为s、下半带宽为r的带状矩阵作选主元LU （Doolittle）分解，将破坏L和U的带状性质；
5：对上半带宽为1、下半带宽为1的三对角带状矩阵，可 以有快速追赶法； 6：对非对角占优的三对角矩阵和拟三对角矩阵作快速追 赶法，可能导致求解不精确甚至求解失败。
u1n
u11 u2n
u22
DU0
,
1
数值分析
A LU
LDU0,
A AT
U
T 0
(
DLT
),
由分解的唯一性，
U
T 0
L.
定理10(对称矩阵的三角分解) 设A为n阶对称矩阵,且A
的顺序主子式均不为零，则A可以唯一分解为
A LDLT ,
其中L为单位下三角阵，D为对角阵.
若A为对称正定矩阵,则
u2 j a1 j l21u1 j
li 2
ai 2
li1u12 u22
a11 a12
a21 a22
an1
an2
a1n 1
a2n ann
l21 ln1
1 ln2
u11 u12
u22
1
u1n u2n unn
数值分析
k 1
比较第k行： akj lkrurj ukj r 1
Ux (14, 10, 72)T x (1, 2, 3)T .
数值分析
L为下三角，U为单位上三角
a11 a12
a21 a22
a n1
an2
a1n l11
a2n ann
l l
21
n1
l22 ln2
1 u12
1
lnn
u1n u12n LU
Ax
b
LUx
b
Ly Ux
a (k ) kk
从矩阵理论来看，相当于左乘矩阵 指标为k的初等下三角阵
1
1
l (k ) k 1k
1
,
li(kk
)
ai(kk ) a(k)
kk
,i
k
1,
,n
l(k) nk
1
数值分析
因此，整个Gauss消元法相当于左乘了一个单位下三角阵
1 l21 L ln1
n
n3/3
GaussJordan elimination
n3/2
Square root
/improved square root
n3/6
追赶法
5n-4
Row/column/ complete Pivoting
Only eliminating elements in a
column below the diagonal one
当需选主元时，PA=LU.
设第r-1步已完成,就有
u11 u12 u1,r1 u1r u1n
l21
u22 u2,r1
u2r
u2n
A lr1,1 lr1,2
ur1,r1 ur1,r
ur
1,n
lr1 lr2 lr,r1
arr
arn
ln1 ln2 ln,r1 anr ann
Cholesky factorization
LU factorization
Jacobi method Gauss-Seidel method
SOR method GCR
In Scientific Computing
↓
Large Linear Systems
Ax=b
as sub-problems/ as intermediate steps
1
l(k) k 1k
1
l(k) nk
l (n) nn1
1
所以有 L s.t. LA U L为单位下三角阵，U为上三角阵 L s.t. A LU
数值分析
Ax b
A LU
LUx b Ux y
三角方 程组， 易于求
解
Ly b Ux y
矩阵的直接三角分解(LU 分解)对解线性方程组有
什么帮助？
A三对 角，弱
对角占 优
数值分析
直接法内容总结
1：用LU分解（Doolittle）求解线性代数方程组等价于顺序
Gauss 消元法：Ax＝b ↔ LUx＝b ↔ Ux＝y，Ly＝b；
2：用选主元LU分解（Doolittle）求解线性代数方程组等价于选
列主元Gauss消元法：Ax＝b ↔ LUx＝Qb ↔ Ux＝y，Ly＝Qb；
北京航空航天大学 数学与系统科学学院
Email: numerical_analysis@hotmail.com Password:beihang 答疑时间：星期四下午2:30－5:30 答疑地点：主216
朱立永
数值分析
第二章 线性方程组的解法
第四讲 矩阵的三角分解法
数值分析
Gaussian elimination
)
2
,
k 1
j 1
称为平方根法，
因为带了开方运算， 因此不常用
2. lij (aij likl jk )/l jj , (i j 1, , n).
k 1
解AX=b的平方根法:…
数值分析
Summary
Cramer rule
(n+1)!
Gauss elimination
n3/3
LU factorizatio