li2?ai2?li1u12u22k?1i?1,i?3,?,n?5?
一般地,设已经给出U的第1行到第k行元素为
?1行元素与L的第1列到第k?1列元素，则U的第kukj?akj??lkiuij,j?k,k?1,?,n?6?L的第k列元素为
lik?aik??lijujkj?1k?1
ukk,i?k?1,k?2,?,n?7?
总结上述讨论,可得用直接三角分解求L,U的计算公式
j?1,2,...n?u1j?a1j,??l?ai1,i?2,3,...ni1?u11?k?1?ukj?akj??lkiuij,j?k,k?1,...,n,k?2,3...,ni?1?k?1?aik??lijujk?j?1l?,i?k?1,...,n?ikukk?
uaula由于在电算时，当kj计算好后,kj就不用了,而ik算好后ik也不再使用,因此,计算好的kj与lik可以存放在A的相应位置,例如
?a11a12a13a14??u11u12?aa22a23a24??l21u2221???A???a31a32a33a34??l31l32???aaaa?41424344??l41l42最后,在存放A的数组中,得到分解矩阵L,U的元素. 例2 例2 将矩阵A进行LU分解
?4215??87210??A???4836???1261120??
解
u13u23u33l43u14?u24??u34??u44?
?4?8A???4??12215??47210??8???836??4??61120??12215?7210??836??61120?
?4?2???1??3215??427210??23???836??18??61120??36?4215??4?2300??2??????1236??1???301120???3?4215??4?2300??2??????1221??1???30420???3?1?2L???1??3所以
15?00??36??1120? 215?300??221??01120? 215?300??221??041?
000??4215??0300?100??,?U???0021?210????041?0001??
下面考虑求解矩阵的LU分解及由此求解方程组Ax?b的运算量.第1步要进行除法n?1次，第2步要进行除法n?2次，乘法与加法均为?n?1???n?2?次；一般地.第k步要进行除法n?k次,乘法与加法?k?1??n?k???k?1??n?k?1?次，所以LU分解共
需加减法
k?1次，乘除法
1????????k?1n?k?k?1n?k?1?n?n3337?n2?n26
3k?1?n?k?1??n?k???k?1??n?k?1???n?k???1n3次.而由此求解
Ax?b还需进行加减法
k?12???k?1???n?k???n?nn2?n2?n3
次;乘除法
k?12???k?1???n?k??1??nn
13213121n?nn?n?n263；次.故由LU分解求解线性方程组的运算量乘除法为3加减法3次.与利用高斯消元法的计算量基本相同. 从直接三角分解公式可看出，当
ukk?0时，计算将中断或者当ukk绝对值很小时,按分解公
式计算可能引起舍入误差的累积.因此，对非奇异矩阵
A可采用与列主元素消元法类似的方法，
将直接三角分解法修改为列主元三角分解法.
?1步分解已完成，这时有
?u11u12?lu22?21???A??uk?1k?1uk?1k?lkk?1akk?????ank?ln1ln2?lnk?1u为了避免用小的数kk作除数,第k步分解需引入
设第ku1n?u2n??????uk?1n??akn?????ann??
si?aik??lijujk,j?1k?1i?k,k?1,?,nsi,skk?i?n?9??10?于是有
ukk?sk,令
lik?i?k?1,?,n
sik?maxsi则交换
A的第k行与ik行元素,然后再进行第k步分解计算,于是有
lik?1,i?k?1,...,nL1P1iA?0??A?1?1
下面用矩阵运算来描述列主元素法为
?11? ?12?
?13?
LkPkiA?k?1??A?k?,kk?1,2,...n?1k其中第
Lk的元素满足lik?1,i?k?1,...n,Pki是初等排列矩阵(由交换单位矩阵的第k行与
ik行得到，从而
Ln?1Pn?1i...L2P2iL1P1iA?0??A?n?1??Un?121简记为L~?14?
?U,其中
L?Ln?1Pn?1i...L2P2iL1P1in?12~1令
?15?
Ln?1?Ln?1
Ln?2?Pn?1iLn?2Pn?1in?1~~n?1
n?2n?1Ln?3?Pn?1iPn?2iLn?3Pn?2iPn?1in?1n?2~ ……
?16?
n?1L1?Pn?1iPn?2i...P2iL1P2i...Pn?2iPn?1in?1n?222n?2~
则Lk是单位下三角矩阵，且
~Ln?1Ln?2...L1Pn?1iPn?2i...P1in?1n?2~~~1
~1记
?Ln?1Pn?1iPn?2i...L2P2iL1P1i?L
n?1n?22P?Pn?1iPn?2i...P1in?1n?21?17?
L?Ln?1Ln?2...L1则
?1~~~?18?
?19?
PA?LU其中P为排列矩阵,L为单位下三角矩阵,U为上三角矩阵,总结以上的讨论有 定理2(列主元素的三角分解定理) 如果A为非奇异矩阵，则存在排列矩阵P,使得
PA?LU其中L为单位下三角阵, U为上三角阵.
在列主元素的三角分解中，L的元素存放在数组A的下三角部分；U的元素存放在A的上三
角部分，而P可通过整型数组P(n)表示. 三 平方根法
应用有限元法解结构力学问题时,最后归结为求解线性方程组,系数矩阵大多具有对称正定性，所谓平方根法就是利用对称正定矩阵的三角分解而得到的求解对称正定方程组的一种有效方法.目前在计算机上广泛应用平方根法解此类方程组.
不难证明,在满足定理 1的条件下,有下面结果.
定理3(矩阵的LDU分解) 设A为n阶矩阵，如果
A的各阶顺序主子式，
Di?0,i?1,2,...n，则A可唯一地分解为
A?LDU其中L为单位下三角阵,U为单位上三角阵,D为对角阵.进一步,如果A是对称矩阵,则
?20?
U?LT,即
A?LDLT
定理4(对称正定矩阵的三 角分解) 如果A为n阶对称正定矩阵,则存在一个实的非奇异下三
角阵L，使得
A?LLT当限定L的对角元素为正时,这种分解是唯一的.
由矩阵乘法,可直接获得L的计算公式
2?lkk???akk??lkj?j?1??k?112?21?
?22?i?k?1,...n
lik?(aik??lijlkj)j?1k?1lkkk?1,2,...n按此方法进行的矩阵分解称为平方根法.由于
2akk??lkj,j?1k?23?
k?1,2,...n所以
2lkj?akk?max?akk?1?k?n于是
lkj??max?akk?max?2k,j1?k?n因此,分解过程中元素kj的数量级不会增长,且对角元素
换,不选主元的平方根法就是一个数值稳定的方法. 例3 例3 用平方根法分解对称正定矩阵
l
lkk恒为正数,于是无需进行行的交
?11??4A???14.252.75?????12.753.5??
l11?a11?4?2l21?解
a21?1???0.5l112 a311??0.5l112
l31?2l22?a22?l21?4.25?0.25?2
a?ll2.75?0.5??0.5?l32?323121??1.5l222
22l33?a33?l31?l32?3.5?0.25?2.25?1
T于是A?LL，其中
00??2L???0.500?????0.51.51??
1n?n?1?由于A为对称矩阵,因此,在电算时只要存储A的下三角部分,其需要存储2个元素,
可用一维数组存放,即
?1?A?n?n?1????a11,a21,...,an1,an2,...ann??2? ?1?1A?n?n?1??i?i?1??jaij?的第2矩阵元素存放在?2个位置,L的元素存放在A的相应位
置上.另外,平方根法的运算量是
开平方 n次；
13321n?n?n23次； 乘除法 6137n?n2?n6次. 加减法 6