所以,
ukj?akj??lkrurjr?1nj?k,k?1,...,nnl?(aik??lkrurj)/ukk同理可推出计算L的第k列的公式: ikr?1i?k,k?1,...,n因此得到如下算法——杜利特(Doolittle)算法:
(1)将矩阵分解为A=LU,对k=1,2,…,n;j=k,k+1,…n; i=k,k+1,…n;
n? ?uki?akj??lkiuyjr?1?n??ljk?(aik??lkyurj)/ukk公式1 r?1??lkk?1??
(2)解Ly=b
公式2:yk?bk??lkryrr?1k?1
k?1,2,...,n(3)解Ux=y
公式3:xk?(yk?r?k?1?unkrrx)/ukkk?n,n?1,...,1对大规模稀疏问题,如果能够通过调整方程及未知量的顺序使得方程组的系数矩阵成带状结构,则对系数矩阵使用通常的LU分解,可以保障单位下三角矩阵L及上三角矩阵U仍为带状结构.
3、直接三角分解法
Gauss消去法还有许多变形,有些变形是为了利用特殊技巧减少误差,把Gauss消去法改写为更紧凑的形式,还有一些变形时根据某类矩阵的特性作一些修正和简化,这些方法可统称为直接三角分解法。
矩阵的三角分解
设A的顺序主子式?k?0(k?1,2,L,n),则可建立线性方程组Ax?b的Gauss消去法与矩阵分解的关系,即矩阵A的LU分解。这个问题前面已经讲的比较详
细了,此处不再赘述。
Doolittle分解法
首先假设A的顺序主子式?k都不为零,则A可作Doolittle分解,即
A?LU,其中L是单位下三角阵,有lii?1,i?j时lij?0;U是上三角阵,i?j时uij?0。仔细写出为
?a11a12?aa22A??21?LL??an1an2LLLLa1n??1??a2n??l211?L??LLO??ann??ln1ln2L??u11u12L??u22L????O??1??u1n??u2n??LU(2.11) L??unn?(k)
在前面逐步推导L和U的元素公式都要借助于有关的aij来表示。现在强调
指出,只要从给定的A通过比较(2.11)式的两边就可能逐步地把L和U构造出来,而不必利用Gauss消去法的中间结果,这种方法称为Gauss消去法的紧凑格式。
根据矩阵的乘法规则,比较(2.11)式的两边可得
min(i,j)aij??k?1likukj,i,j?1,2,L,n (2.12)
先算出U的第1行,在(2.12)令i?1,由l11?1,l1j?0(1?j?n)可得
u1j?a1j,j?1,2,L,n
接着在(2.12)中令j?1,由ai1?li1?u11,从而算出L的第1列
li1?ai1/u11?ai1/a11,i?2,L,n
若U的第1至k?1行和L的第1至k?1列已经算出,则由(2.12)可计算出
U的第k行元素
ukj?akj??lkrurj,j?k,k?1,L,n (2.13)
r?1k?1接着再计算出L的第k列元素
lik?aik??lirurkr?1k?1ukk,j?k?1,L,n (2.14)
解方程组Ax?b就化为求解LUx?b,分两步求解。第一步解Ly?b,其中L为下三角阵,只要用逐次向前代入的方法;第二步解Ux?y,其中U为上三角阵,用逐次向后代入方法即可解除x。
例 用Doolittle方法求解:
?6??2?1???1?1??x1??6??????410??x2???1? ??????14?1x35?????0?13??x4???5?21解:第1步算出L和U:
?1??11?3L??11?5?6?11???610第2步解Ly?b得:
??62??10????3?,U??1??????9?1??37??1233710?1??1?3?9? ??10?191??74?23191??y??6,3,,??
574??第3步解Ux?y得:
x?(1,?1,1,?1)T
T
二、 迭代法
迭代法具有的特点是速度快。与非线性方程的迭代方法一样,需要我们构造一个等价的方程,从而构造一个收敛序列,序列的极限值就是方程组的根。 对方程组:Ax?b做等价变换:x?Gx?g 如:令A?M?N,则
Ax?b?(M?N)x?b?Mx?b?Nx?x?M?1Nx?M?1b
则,我们可以构造序列x(k?1)?G x(k)?g
若x(k)?x*?x*?G x*?g?Ax*?b 则,我们可以构造序列x(k?1)?G x(k)?g 若x(k)?x*?x*?G x*?g?Ax*?b 同时:x(k?1)?x*?Gx(k)?Gx*?G(x(k)?x*)
?L?Gk?1(x(0)?x*)
所以,序列收敛?Gk?0 与初值的选取无关 1. Jacobi迭代
?a11x1?L?a1nxn?b1?M ??ax?L?ax?bnnnn?n11?1?x??1a(a12x2?L?a1nxn?b1)11??1?x?(a21x1?a23x3?L?a1nxn?b2)?2a22 ? ??M??1?x??na(an1x1?L?an n?1xn?1?bn)nn??(k?1)?x1??(k?1)?x2? ????(k?1)?xn????1(a12x2(k)?L?a1nxn(k)?b1)a11?1(a21x1(k)?a23x3(k)?L?a1nxn(k)?b2)a22 ?1(an1x1(k)?L?an n?1xn?1(k)?bn)annM?格式很简单:
xi(k?1)?1i?1?(?aijxj(k)?aiij?1j?i?1?axijn(k)j?bi)
迭代矩阵
记A?D?L?U
?a11?D??O?0?0??? ann???0???a210L??MO????a?n1OL0?ann?10??0?a12L??0O??? U??O??????00???a1n??M?? ?0?an?1n?0??易知,Jacobi迭代有
(D?L?U)x?b Dx?(L?U)x?b x?D?1(L?U)x?D?1b
? B0?D?1(L?U)?I?D?1A , f0?D?1b,B0是简单迭代的迭代矩阵。
看上述公式和过程很抽象,我们来举个简单例子:
?a11x1?a12x2?a13x3?b1??a21x1?a22x2?a23x3?b2 (aii?0) ?ax?ax?ax?b?3113223333?1x?(b1?a12x2?a13x3)?1a11??1得?x2?(b2?a21x1?a23x3)
a22??1(b3?a31x1?a32x2)?x3?a33?(0)(0)T,x3)作初始近似,代入, 得上述变换的方程后,我们任取一向量x(0)?(x1(0),x2(1)(1)T(m)(m)T,x3),L,x(m)?(x1(m),x2,x3) x(1)?(x1(1),x2 假设上述方程的准确解是a?(a1,a2,a3)T
那么如果limxm?a,?我们就说上述迭代是收敛的。
m??2. Gauss-Seidel迭代
在Gauss-Seidel迭代中,使用最新计算出的分量值