线性方程矩阵三角分解法

1、计算物理计算物理矩阵的直接三角分解法矩阵的直接三角分解法列主元三角分解法列主元三角分解法11计算物理计算物理矩阵三角分解法矩阵三角分解法Gauss消元法：增广矩阵（A|b）经过若干次行初等变换，使之化 为（U|g）形式，其中矩阵U为上三角矩阵由线性代数理论可知，对一个矩阵进行一次初等行变换，相当于该矩阵左乘一个相应的初等矩阵，因此Gauss消元过程可表示为)|()|(.121gUbALLLLkk其中 为单位下三角矩阵，即),.2, 1(kiLi若记 （仍为下三角矩阵），则 1121).(LLLLLkkLUA UALLLLkk121.Gauss过程隐含了系数矩阵的一个矩阵三角分解计算物理计算物理

2、三角分解定义三角分解定义设A是n阶方阵（n=2），则称A=LU为矩阵A的一个三角分解，其中L是下三角矩阵，U是上三角矩阵Doolittle分解L是单位下三角矩阵，U是上三角矩阵Crout分解L是下三角矩阵，U是单位上三角矩阵定理定理若n阶矩阵A的k阶顺序主子式不为零，则A有唯一的Doolittle分解与Crout分解)1,.2, 1(nk计算物理计算物理nnnnnnnnnnnnnnnnnnnnnnnuCBuBCCLCLUBBUULULuCUBLuCUBLaAAAAULLUA) 1() 1() 1() 1() 1() 1() 1() 1() 1() 1() 1() 1(21) 1(110) 1(

3、nnnnnnnnnuuCCBBUULL) 1() 1() 1() 1(可逆0000)1()1()1()1(nnnnUULL1)1()1()1(1)1(nnnnUULL11)1()1()1(1)1(nnnnnIUULL计算物理计算物理直接三角分解法直接三角分解法设矩阵A存在Doolittle分解A=LU，则方程组bxLUbAx)(yUxbLynnnnbbbyyylll21212121111nnnnnnyyyxxxuuuuuu212122211211计算物理计算物理nnnnbbbyyylll21212121111nnnnnnyyyxxxuuuuuu2121222112111111kmmkmkkyl

4、byby)(11nkmmkmkkkknnnnxuyuxuyxnk,.3 , 21,.2, 1nnk计算物理计算物理为得到方程组的解，需先求出矩阵A的Doolittle分解LUA nnnnnnnnnnnnuuuuuulllaaaaaaaaa222112112121212222111211111由矩阵乘法得111111111111ualulaauuaiiiijjjj),.2,1(nj U的第一行L的第一列),.3,2(ni 计算物理计算物理),.1,(11nkkjuulakmkjmjkmkj),.2,1(11nkkiululakmkkikmkimikU的第k行、L的第k列),.1,(11nkkju

5、laukmmjkmkjkj),.2,1()(111nkkiulaulkmmkimikkkik),.3,2(nk 计算物理计算物理为了节约存储空间，方便编程计算，将每一步中的结算结果u、l值存放在对应的a所占的空间位置nnnnnnnnnnnnulluuluuuaaaaaaaaaA2122221112112122221112111111kmmkmkkylbyby)(11nkmmkmkkkknnnnxuyuxuyx计算物理计算物理直接三角分解法解题步骤先求U矩阵的行元素，再求L矩阵的列元素一行一列交叉进行计算，将计算值存放到矩阵A的对应位置得到完整的U、L矩阵系数矩阵A计算列向量y得到方程组解x计算

6、物理计算物理列主元三角分解法列主元三角分解法目标：与Gauss列主元消元法的原理相同，为了消除因除数过小而产生的误差扩散),.1,(11nkkjulaukmmjkmkjkj),.2,1()(111nkkiulaulkmmkimikkkik第K步),.2,1(11nkkiulaskmmkimikik11kmmkkmkkkkkkulaus),.2,1(nkkisslkkikik计算物理计算物理为了消除误差扩散，应当满足ikkkss但分解公式中不一定满足该条件，此时选行号 ，使kiiknikkissk max若 ，则交换矩阵A中第k行和第 行的元素 kikki此时一定满足ikkksskks),.1,

7、(nkkisik（ 是 中的主元）计算物理计算物理直接三角分解法解题步骤求U矩阵的第一行元素与L矩阵的第一列元素一行一列交叉进行计算，将计算值存放到矩阵A的对应位置得到完整的U、L矩阵增广矩阵A|b计算列向量y得到方程组解x列主元三角分解法解题步骤第一步：求出第一列的主元 ，交换行11ia第k步：求出主元 ，交换行kiks求U矩阵的第k行元素与L矩阵的第k列元素计算物理计算物理例：用列主元三角分解法解下列方程组553618454622321321321xxxxxxxxxMatlab程序计算物理计算物理a=2,1,2;4,5,4;6,-3,5;b=6;18;5;n=3;A=a,b;u=;l=;s

8、=;y=;x=;p=A(1,1);q=1;for i=1:n if p=A(i,1) p=p; q=q; else p=A(i,1); q=i; end;end;for j=1:n+1 t=A(1,j); A(1,j)=A(q,j); A(q,j)=t;end;for j=1:n u(1,j)=A(1,j); A(1,j)=u(1,j);end;for i=2:n l(i,1)=A(i,1)/u(1,1); A(i,1)=l(i,1);end;for k=2:n for i=k:n d=0; for m=1:k-1 d=d+A(i,m)*A(m,k); end; s(i,k)=A(i,k)-d

9、;end; p=s(k,k); q=k; for i=k:n if p=s(i,k) p=p; q=q; else p=s(i,k); q=i; end; end; for j=1:n+1 t=A(k,j); A(k,j)=A(q,j); A(q,j)=t; end；计算物理计算物理for j=k:n d=0; for m=1:k-1 d=d+A(k,m)*A(m,j); end; u(k,j)=A(k,j)-d; A(k,j)=u(k,j); end; if k=n-1 for i=k+1:n d=0; for m=1:k-1 d=d+A(i,m)*A(m,k); end; l(i,k)=1/u(k,k)*(A(i,k)-d); A(i,k)=l(i,k); end; end;end;y(1,1)=A(1,n+1);for k=2:nd=0; for m=1:k-1 d=d+A(k,m)*y(m,1); end; y(k,1)