计算方法线性方程组的直接解法

第六章线性方程组的直接法2/109§2.0引言§2.1Gauss消去法§2.2矩阵分解在解线性方程组中的作用补充向量的范数与矩阵的范数§2.3直接法的误差分析§2.4线性方程组迭代解法

3/109§2.0引言在自然科学和工程技术中很多问题的解决常常归结为解线性代数方程组。例如电学中的网络问题，用最小二乘法求实验数据的曲线拟合问题，解非线性方程组问题，用差分法或者有限元方法解常微分方程、偏微分方程边值问题等都导致求解线性代数方程组。而这些方程组的系数矩阵大致分为两种，一种是低阶稠密矩阵(例如，阶数大约为≤150)，另一种是大型稀疏矩阵(即矩阵阶数高且零元素较多)。

4/109§2.0引言设有线性方程组Ax=b，其中为非奇异阵，，关于线性方程组的数值解法一般有两类：直接法与迭代法。5/1091.直接法就是经过有限步算术运算，可求得方程组精确解的方法(若计算过程中没有舍入误差)。但实际计算中由于舍入误差的存在和影响，这种方法也只能求得线性方程组的近似解。本章将阐述这类算法中最基本的高斯消去法及其某些变形。这类方法是解低阶稠密矩阵方程组的有效方法，近十几年来直接法在求解具有较大型稀疏矩阵方程组方面取得了较大进展。6/1092.迭代法迭代法是用某种极限过程去逐步逼近线性方程组精确解的方法。迭代法具有需要计算机的存贮单元较少、程序设计简单、原始系数矩阵在计算过程中始终不变等优点，但存在收敛性及收敛速度问题。迭代法是解大型稀疏矩阵方程组(尤其是由微分方程离散后得到的大型方程组)的重要方法。7/109§2.1Gauss消去法高斯（Gauss）消去法是解线性方程组最常用的方法之一。它的基本思想是通过逐步消元（行的初等变换），把方程组化为系数矩阵为三角形矩阵的同解方程组，然后用回代法解此三角形方程组（简单形式）得原方程组的解。

8/109§2.1Gauss消去法例：下面我们来讨论一般的解n阶方程组的高斯消去法。9/1091.消元将原方程组记为A(1)x=b(1)，其中A(1)=(aij(1))=(aij)，b(1)=b(1)第一次消元。其中注：若a11(1)=0，则在第1列中至少有一个元素不为0，可交换行后再消元10/109(2)第k次消元。其中注：为减少计算量，令，则11/109(3)当k=n

–1时得完成第n

–1次消元后得到与原方程组等价的三角形方程组A(n)x=b(n)注：当det(A)≠0时，显然有aii(i)≠0，(i=1，…，n)，称为主元素。12/1092.回代求解三角形方程组A(n)x=b(n)，得到求解公式注：求解过程称为回代过程。13/1093.Gauss消去法的计算量以乘除法的次数为主(1)

消元过程：第k步时(n–k)+(n–k)(n–k+1)=(n–k)(n–k+2)共有14/1093.Gauss消去法的计算量(2)

回代过程：求xi中，乘n–i次，除1次，共n–i+1次（i=1，…，n–1）共有15/1093.Gauss消去法的计算量(3)总次数为注：当n=20时约为2670次，比克莱姆法则9.71020次大大减少。16/1094.说明(1)若消元过程中出现akk(k)=0，则无法继续(2)若akk(k)≠0，但较小，则小主元做除数将产生大误差(3)每次消元都选择绝对值最大者作主元，称为高斯主元消去法(4)通常第k步选择——第k列主对角元以下元素绝对值最大者作主元（该行与第k行对调），称为列主元消去法。17/109例1用舍入三位有效数字求解线性方程组（准确解是x1=10.0，x2=1.00）解：1)

不选主元的Gauss消去法计算结果：

x1=-10.0，x2=1.01，此解无效；

2)

按列选主元的Gauss消去法计算结果：

x1=10.0，x2=1.00.18/109(5)行标度化当线性方程组的系数矩阵的元素大小相差很大时，则可能引进因丢失有效数字而产生的误差，并且舍入误差的影响严重，即使用Gauss主元消去法得到的解也不可靠．为避免这一问题，可将系数矩阵的行元素按比例增减以使元素的变化减小．如用每行元素的最大模除该行各元素，使它们的模都不大于1，这叫行标度化，其目的是要找到真正的主元．19/109(5)行标度化如用每行元素的最大模除该行各元素，使它们的模都不大于1，这叫行标度化，其目的是要找到真正的主元．消元过程仍是对原方程组进行的，只不过在Gauss列主元消去法的算法中，按列选主元ck时，应修改为其中20/1095.算法设计输入A，n，epsFor(k=1ton-1)选主元A(I0，k)——确定行号p=A(I0，k)若|p|<=epsText=1，退出循环若I0kT换行（I0k）消元若|A(n，n)|<=epsText=1F回代若ext=1T无解F输出解I0=kfor(i=k+1ton)若|A(i，k)|>|A(I0，k)|TI0=ifor(i=k+1ton)A(i，k)=A(i，k)/A(k，k)for(j=k+1ton+1)A(i,j)=A(i,j)-A(i,k)*A(k,j)A(n，n+1)=A(n，n+1)/A(n，n)for(k=n–1to1step-1)w=0for(j=k+1ton)w=w+A(k，j)*A(j，n+1)A(k，n+1)=A(k，n+1)–wA(k，n+1)=A(k，n+1)/A(k，k)for(j=kton)n=A(k，j),A(k，j)=A(I0，j),A(I0，j)=n21/1095.算法设计(1)

高斯消去法a=[5-1-1-1-4;-110-1-112;-1-15-18;-1-1-11034];x=[0;0;0;0]n=4fork=1:n-1fori=k+1:na(i,:)=a(i,:)-a(k,:)*a(i,k)/a(k,k)endendforj=n:-1:1

x(j)=(a(j,n+1)-a(j,j+1:n)*x(j+1:n))/a(j,j)end22/109(2)列主元素消去法a=[-0.040.040.123;0.56-1.560.321;-0.241.24-0.280];x=[0;0;0];n=3;fork=1:n-1[c,i]=max(abs(a(k:n,k)));

q=i+k-1ifq~=km=a(q,:);a(q,:)=a(k,:);a(k,:)=mendfori=k+1:na(i,:)=a(i,:)-a(k,:)*a(i,k)/a(k,k)endendforj=n:-1:1

x(j)=(a(j,n+1)-a(j,j+1:n)*x(j+1:n))/a(j,j)end23/109§2.2矩阵分解在解线性方程组中的作用1.原理若A=LR，则Ax=b

LRx=b

若其中L、R为三角形矩阵，则方程组易解定义：(1)称为单位下三角阵(2)设L为单位下三角阵，R为上三角阵，若A=LR，则称A可三角（LR）分解，并称A=LR为A的三角分解（杜利特尔分解）24/109定理：(1)A=(aij)nn有唯一LR分解

A的顺序主子式k≠0，k=1，2，...，n

–1(2)若A=(aij)nn可逆，则存在置换阵P，使PA的n个顺序主子式全不等于0注：由Ax=b

PAx=Pb知，经适当行交换后，A总可三角分解25/1092.直接三角分解法设A的各阶顺序主子式不为0，且A=LR，即(1)主对角线（含）上边：当列标m

>

行标k时，lkm=0

j=k，k+1，...，n；k=1，2，...，n(2)主对角线下边：当行标m

>

列标k时，rkm=0

i=k+1，k+2，...，n；k=1，2，...，n–126/109设A的各阶顺序主子式不为0，且A=LR，即(1)主对角线（含）上边：当列标m

>

行标k时，lkm=0

j=k,k+1,...,n；k=1,2,...,n(2)主对角线下边：当行标m>列标k时，rkm=0

i=k+1,k+2,...,n；k=1,2,...,n–1欲求lik和rkj27/109(1)主对角线（含）上边：当列标m

>

行标k时，lkm=0

j=k,k+1,...,n；k=1,2,...,n(2)主对角线下边：当行标m>列标k时，rkm=0

i=k+1,k+2,...,n；k=1,2,...,n–1欲求lik和rkj设k=1，a1j=l11r1j

r1j=a1j

j=1，2，...，nai1=li1r11

li1=ai1/r11

i=2，3，...，n28/109一般地，最后：r11r12...r1n第1步l21r22...r2n第2步l31l32......ln1ln2rnn第n步29/1093.求解方程组下三角Ly=b：上三角Rx=y：紧凑格式：30/109a=[2100-3;-3-4-1213;123-4;4149-13];b=[10;5;-2;7];y=[0000]’;x=[0000]'n=4;forr=1:na(r,r:n)=a(r,r:n)-a(r,1:r-1)*a(1:r-1,r:n)a(r+1:n,r)=(a(r+1:n,r)-a(r+1:n,1:r-1)*a(1:r-1,r))/a(r,r)endR=triu(a,0)L=eye(n)+tril(a,-1)forr=1:ny(r)=b(r)-L(r,1:r-1)*y(1:r-1)endforj=n:-1:1x(j)=(y(j)-R(j,j+1:n)*x(j+1:n))/R(j,j)end31/109紧凑格式a=[2100-310;-3-4-12135;123-4-2;4149-137];x=[0000]'n=4;forr=1:na(r,r:n+1)=a(r,r:n+1)-a(r,1:r-1)*a(1:r-1,r:n+1)a(r+1:n,r)=(a(r+1:n,r)-a(r+1:n,1:r-1)*a(1:r-1,r))/a(r,r)endR=triu(a,0)forj=n:-1:1x(j)=(a(j,n+1)-R(j,j+1:n)*x(j+1:n))/R(j,j)end32/1094.选主元的三角分解法5.平方根法(1)定理：若A正定，则A可唯一分解为A=LDLT。其中L为单位下三角，D为元素全为正数的对角阵。(2)设A=LDLT则有：33/109(3)为了避免重复计算，引进中间量tik=likdk，将上式改写为：34/109(4)Ax=b

LDLTx=b

Ly=b，Dz=y，LTx=z求解公式：上述方法称为“改进的平方根法”注：不能选主元作行交换——破坏对称性，必是数值稳定的35/109a=[5-41;-46-4;1-46];b=[2;-1;-1]d=[000]';t=[000;000;000];L=[000;000;000];n=3fork=1:n

d(k)=a(k,k)-t(k,1:k-1)*L(k,1:k-1)'

t(k+1:n,k)=a(k+1:n,k)-t(k+1:n,1:k-1)*L(k,1:k-1)'L(k+1:n,k)=t(k+1:n,k)/d(k)endx=[000]';y=[000]';z=[000]';fori=1:ny(i)=b(i)-L(i,1:i-1)*y(1:i-1)endfori=1:nz(i)=y(i)/d(i)endfori=n:-1:1x(i)=z(i)-L(i+1:n,i)'*x(i+1:n)end36/1096.追赶法在实际问题中，经常遇到以下形式的方程组这种方程组的系数矩阵称为三对角矩阵。以下针对该方程组的特点提供一种简便有效的算法——

追赶法。37/109设直接三角分解为：A=LR其中容易看出pi=ci（1，2，…，n

–1），计算L和R的其余元素的公式为38/109容易看出pi=ci（1，2，…，n

–1），计算L和R的其余元素的公式为计算过程：r1→l2→r2→...→ln→rn39/109有了A的LR分解后，线性方程组Ax=d等价于两个简单的方程组：Ly=d，Rx=y于是，计算y的公式是y1=d1，yi=di

–

liyi–1

（i=2，3，…，n）计算过程：r1→y1→l2→r2→y2→...→ln→rn→yn

（追）计算x的公式是xn=yn/rn，xi=(yi–cixi+1)/ri

（i=n–1，n–2，…，1）(赶）40/109定理：若A为对角占优（diagonallydominant）的三对角阵，且满足|b1|>|c1|>0，|bi|≥|ai|+|ci|，aici≠0，i=2，3，…，n-1|bn|>|an|>0.则追赶法可解以A为系数矩阵的方程组注：1.如果A是严格对角占优阵，则不要求三对角线上的所有元素非零。2.根据不等式可知：分解过程中，矩阵元素不会过分增大，算法保证稳定。3.运算量为O(6n)。41/1097.QR算法42/109补充向量的范数与矩阵的范数在线性方程组的数值解法中，经常需要分析解向量的误差，需要比较误差向量的“大小”或“长度”。那么怎样定义向量的长度呢？我们在初等数学里知道，定义向量的长度，实际上就是对每一个向量按一定的法则规定一个非负实数与之对应，这一思想推广到维线性空间里，就是向量的范数或模。用Rn表示n维实向量空间，用Cn表示n维复向量空间，首先将向量长度概念推广到Rn(或Cn)中。43/1091.向量的范数向量的范数可以看作是描述向量“大小”的一种度量.范数的最简单的例子，是绝对值函数：有三个熟知的性质：(1)x

0

x

>0

x

=0当且仅当x=0(2)ax=

a

x

a为常数(3)

x+y

≤

x

+

y

44/1091.向量的范数范数的另一个简单例子是三维欧氏空间的长度设x=(x1,x2,x3)，则x的欧氏范数定义为：欧氏范数也满足三个条件：

x，y

R3，a为常数(1)

x

≥0，且x=0

x=0(2)

ax

=

a

x

(3)

x+y

≤

x

+

y

前两个条件显然，第三个条件在几何上解释为三角形一边的长度不大于其它两边长度之和。因此，称为三角不等式。45/109向量范数的一般概念：定义1：设V是数域F上的向量空间，对V中任一向量α，都有唯一实数α

与之对应，满足如下三个条件：1)正定性：α≥0，且α=0

α=02)齐次性：kα=|k|α

，这里k

F3)三角不等式：α+

α

+

则称α为α的范数。定义了范数的向量空间称为赋范向量空间.简单性质：(1)x

0——单位向量(2)||x||=||–x||(3)|||x||–||y|||||x–y||——当x

y时，||x||||y||46/109Cn上的常见范数有：1)1-范数

2)2-范数称为欧氏范数3)-范数不难验证，上述三种范数都满足定义的条件。注：上述形式的统一：

1

p

47/109定理5：定义在Cn上的向量范数||x||是变量x分量的一致连续函数。（f(x)=||x||）定理6：在Cn上定义的任何两个范数都是等价的。即存在正数k1与k2(k1≥k2>0)，对一切xCn，不等式k1||x||b

||x||a

k2||x||b成立。对常用范数，容易验证下列不等式：48/109例1设A为正定矩阵，xCn，令，则||x||A为向量范数，称为加权范数证明：(1)显然xA≥0，且xA=0

x=0(2)(3)因为A正定，所以可逆P，使A=PTP

所以||x+y||A=||P(x+y)||2=||Px+Py||2||Px||+||Py||2=||x||A+||y||A综上，||x+y||A为范数。注：当0<p<1时，不是范数49/109有了范数的概念，就可以讨论向量序列的收敛性问题。定义2：设给定Cn中的向量序列{xk}，即x0，x1，…，xk，…其中若对任何i(i=1,2,…,n)都有则向量称为向量序列{xk}的极限，或者说向量序列{xk}依坐标收敛于向量x*，记为50/109定义2：设给定Cn中的向量序列{xk}，即x0，x1，…，xk，…其中若对任何i(i=1,2,…,n)都有则向量称为向量序列{xk}的极限，或者说向量序列{xk}依坐标收敛于向量x*，记为定理7：向量序列{xk}依坐标收敛于x*的充要条件是如果一个向量序列{xk}与向量x*满足上式，就说向量序列{xk}依范数收敛于x*，于是便得：向量序列依范数收敛与依坐标收敛是等价的。51/1092.矩阵的范数下面我们给出矩阵范数的一般定义。为简单起见，仅讨论实数域的情形。定义3

若矩阵ARnn的某个非负的实值函数N(A)=||A||，满足条件1)正定性：A≥0，且A=0

A=02)齐次性：kA=|k|||A||

kR3)三角不等式：||A+B||

||A||+||B||4)相容性：||AB||||A||||B||则称N(A)是Rnn上的一个矩阵范数(或模)。52/1092.矩阵的范数定义3

若矩阵ARnn的某个非负的实值函数N(A)=||A||，满足条件1)正定性：A≥0，且A=0

A=02)齐次性：kA=|k|||A||

kR3)三角不等式：||A+B||

||A||+||B||4)相容性：||AB||||A||||B||则称N(A)是Rnn上的一个矩阵范数(或模)。如由Rnn上2-范数可以得到Rnn中矩阵的一种范数称为A的Frobenius范数。||A||F显然满足正定性、齐次性及三角不等式。53/109由于在大多数与估计有关的问题中，矩阵和向量会同时参与讨论，所以希望引进一种矩阵的范数，它是和向量范数相联系而且和向量范数相容的，即5)||Ax||||A||||x||对任何向量xRn及矩阵ARnn都成立。注：满足相容性条件的范数很多，其中最重要的是算子范数54/109再引进一种矩阵的范数。定义4

设xRn及矩阵ARnn，||x||v为一种向量范数(如v=1，2或∞)，相应地定义一个矩阵的非负函数可验证||A||v满足定义3，所以||A||v是Rnn上矩阵的一个范数，||A||v还满足相容性条件5)，称为A的算子范数，或诱导范数。注：(1)可以证明：(2)||Ax||是x的连续函数，从而存在x*使||x*||=1且||Ax*||=||A||(3)若||.||为算子范数，I为单位矩阵，则||I||=155/109定理8：设n阶方阵A=(aij)nn，则与||x||1相容的矩阵范数是与||x||2相容的矩阵范数是其中1为矩阵ATA的最大特征值。与||x||相容的矩阵范数是上述范数分别称为1-范数，2-范数和-范数又称为列模、谱模和行模56/109注：A的范数||A||与A的特征值间的关系设为矩阵A的任一特征值，e为相应的特征向量，则e=Ae，因为||||e||=||Ae||||A||||e||，所以||||A||从而得定理9：矩阵A的任一特征值的绝对值不超过A的范数||A||。定义6：矩阵A的诸特征值的最大绝对值称为A的谱半径，记为：由定理9可知：57/109例(p138)：设A=(aij)nn，||A||为其算子范数，则||A||<1

I–A可逆，且58/109§2.3直接法的误差分析1.误差向量与残差向量设x为线性方程组Ax=b的准确解，用数值算法得到的近似解为x*，则称向量e=x*-x为误差向量。显然，||e||可以表示近似解x*的准确程度：||e||越小，近似解x*越准确。但由于x不知道，所以无法计算e。变通的办法是考虑残差向量：r=b–Ax*的范数的大小。但是，在某些情况下，尽管||r||很小，||e||也可能很大。59/109例线性方程组的准确解是(1，1)T，若用某种方法得到计算解(2.000，0.500)，则残差向量其-范数||r||=0.0015，而误差向量其-范数||e||=1，它是残差向量的667倍。60/109能否建立残差向量与误差向量之间的关系？由于r=b–Ax*=Ax–Ax*=A(x–x*)=Ae，故有e=A–1r，||e||||A–1||||r||另一方面，由b=Ax得||b||||A||||x||，因此其中，是近似解x*的相对误差，而表示相对残差上式表明近似解x*的相对误差大约是相对残差的||A||||A–1||倍。61/1092.条件数与病态方程组定义5：设A为n阶非奇矩阵，称数||A||||A–1||为矩阵A的条件数，记为cond(A)。显然，当A的条件数cond(A)较小时，只要残差向量很小，近似解x*的相对误差就很小，即近似解的精度较高.

但若cond(A)较大时，则即使相对残差较小，近似解x*的相对误差也可能很大，近似解的精度就很低。62/1092.条件数与病态方程组上例中线性方程组系数矩阵A的条件数：||A||=3，||A-1||

1000cond(A)3000

很大！这就解释了为什么残差向量小未必保证近似解的误差小。定义6：当cond(A)相对很大时，称方程组Ax=b为病态的，否则称为良态的。63/1093.病态方程组的判断理论上虽然可以用矩阵的条件数来度量线性方程组求解问题的病态程度.但是，条件数大到怎样才算大，并没有客观的适用的标准，且当方阵的阶n较大时，A-1的计算也非易事，所以一般地说，遇到下列几种情况之一就应考虑线性方程组可能是病态的：

(1)系数矩阵的某两行(列)几乎近似相关；

(2)系数矩阵的行列式的值很小；

(3)用主元消去法解线性方程组时出现小主元；

(4)近似解x*已使残差向量r=b-Ax*的范数很小，但该近似解仍不符合问题的要求．64/1094.解的精度估计与改善我们知道，影响计算精度的一个重要因素是舍入误差的传播．如前所述，按计算过程逐步地分析舍入误差的传播是极为困难又不实用的．一种常用的办法是向后误差分析方法，它将计算过程中舍入误差的影响归结为原始数据的小扰动，也就是说，将计算所得的近似解x*看成是某个带有小扰动的方程组(A+A)x=b+b的准确解．65/109定理：设x*是线性方程组(A+A)x=b+b的解，x是线性方程组Ax=b的解，且A的扰动A满足||A-1||||A||<1，则有上式表明了近似解x*的相对误差与A及b的相对误差之间的关系，其中A的条件数起着重要的作用．66/109估计式只是理论性的，实际应用有很大的困难．因此，—般采用下述方法估计近似解的精确程度．对于线性方程组Ax=b，任取一个非零向量y，计算出d=Ay；然后采用与解Ax=b同一算法解线性方程组Ay=d，得到它的近似解y*，y*的相对误差为由于求解Ay=d与求解Ax=b的算法相同，且系数矩阵同为A，所以可以认为Ax=b的近似解x*的相对误差大约也是当然，对病态方程组来说，这种估计方法是不大可靠的．67/109

病态方程组的求解问题是计算方法研究的一个重要课题．这里介绍一下改善近似解的精度的方法：

1)采用双倍字长进行运算，特别是在做形如∑aibi的计算时，每一乘积项均取双倍字长的数，并与其他项累加后再舍入成单字长，这样能使结果的精度提高，同时也比全都用双倍字长进行运算节省内存和机时．

68/1092)采用迭代改善的办法．设已得到Ax＝b的近似解x(0)，采用双倍字长计算Ax(0)后再舍入成单字长，并计算残差向量r(0)=b–Ax(0)，记x=x–x(0)，则x应满足方程组Ax＝r(0)．（Ax＝A(

x–x(0))＝Ax–Ax(0)=b–Ax(0)=r(0)）采用原单字长运算的算法求出x的近似解x(0)(在A=LR或A=QR分解已完成下，只需作顺代和回代求解)

从而得到Ax＝b的改善了的近似解x(1)=x(0)+x(0)，若||x(0)||<，则x(1)就是满足精度（）要求的解。否则再重复上述做法。69/109§2.4线性方程组迭代解法1.基本思想若方程组Ax=b可写成等价形式x=Bx+g，(2.4-1)则给定一个初始向量x(0)，可以得到迭代公式x(k+1)=Bx(k)+g，k=0，1，…(2.4-2)若上式确定的向量序列{x(k)}收敛于x，则x显然是(2.4-1)的解，从而为方程组Ax=b的解。形如(2.4-2)的逐次逼近的方法称为简单迭代法，B称为该迭代法的迭代矩阵。70/1092.迭代法的收敛性

定理1：对任意初始向量x(0)及常向量g，迭代法(2.4-2)收敛的充分必要条件是迭代矩阵B的谱半径(B)<1。这一结论在理论上是颇为重要的，但实际用起来不甚方便。当n较大时，谱半径的计算并非易事，但由于(B)||B||，所以只要B的某种范数小于1，迭代法(2.4-2)必收敛：

定理2：若迭代矩阵Ｂ的某种范数||B||<1，则(2.4-2)确定的迭代法对任意初值x(0)均收敛于方程组x=Bx+g的唯一解x。71/109

定理3：若||B||=q<1，则(2.4-2)确定的向量序列满足：(2.4-3)(2.4-4)其中x是方程组(2.4-1)的准确解。证明：(1)因为x(k)–x=B(x(k–1)–x)=B(x(k–1)–x(k)+x(k)–x)，所以||x(k)–x||=||B[(x(k–1)–x(k))+(x(k)–x)]||

||B||||(x(k–1)–x(k))+(x(k)–x)||

q(||x(k–1)–x(k)||+||x(k)–x||)，从而得(2.4-3)。72/109(2)因为x(k)–x(k–1)=B(x(k–1)+g)–B(x(k–2)+g)=B(x(k–1)–x(k–2))=…=B

k–1(x(1)–x(0))，所以有||x(k)–x(k–1)||=||B

k–1(x(1)–x(0))||

||B

k–1||||x(1)–x(0)||

||B||

k–1||x(1)–x(0)||代入(2.4-3)即得(2.4-4)(2.4-3)与(2.4-4)给出||x(k)–x||的估计分别称为后验估计与先验估计。73/109先验估计可以在求出x(1)后就可以估计出在精度要求为下大约需要的迭代次数k。从证明过程看出，先验估计是由后验估计放大得到，所以较粗。从(2.4-4)可以知道，q=||B||越小，{x(k)}收敛越快。x(1)–x(0)=Bx(0)+g–x(0)=g–(I–B)x(0)的范数越小，或x(0)满足方程组Ax=b的程度越好，则x(k)与x就越接近。74/109后验估计用计算过程中前后相继两个迭代向量之差的范数来估计后一个近似值的误差，较为可靠，并给出迭代终止准则，例如在q

0.5的条件下，有(1)绝对误差准则当||x(k)–x(k–1)||<时终止(2)相对误差准则当||x(k)–x(k–1)||/||x(k)||<时终止注意当q=||B||接近于1时，此准则可能失效。75/1093.简单迭代法、Jacobi迭代与Seidel迭代若将方阵A分裂为A=M

–

N其中M可逆，则线性方程组Ax=b可写成等价形式Mx=Nx+b，或x=Bx+g其中B=M–1N，g=M–1b称x(k+1)=Bx(k)+g，k=0，1，…

为简单迭代法输入矩阵A，向量b，初值x0，eps计算B，g计算x1=B*x0+gwhilenorm(x1-x0,inf)>epsx0=x1k=k+1x1=B*x0+g输出k，x176/109若有迭代公式（分量式）(2.4-2')由于在计算新分量xi(k+1)时，x1(k+1)，x2(k+1)，…，xi-1(k+1)都已经算出，故可将上式改为(2.4-7)称(2.4-7)为Seidel迭代法（迭代公式）。77/109将式中x1(k+1)，x2(k+1)，…，xn-1(k+1)项移到左边：即有78/109Seidel迭代法的矩阵形式为：x(k+1)=BSx(k)+h其中，BS称为Seidel迭代的迭代矩阵。79/109例1设线性方程组问p，q满足怎样条件时，对任意初始向量，简单迭代和Seidel迭代收敛？解：因为简单迭代法的迭代矩阵，而B的特征多项式：(p–)2=q2，特征值为

=p

iq，所以当|

|=p2+q2<1时，简单迭代法对任意初始向量x(0)收敛。80/109Seidel迭代公式是即其迭代矩阵是BS的特征多项式是从而根据E.I.Jury准则，若p2<1，q2<(1+p)2，则有(BS)<1，即有Seidel迭代对任意初始向量x(0)收敛。81/109例2用简单迭代法解线性方程组问用后验估计与先验估计迭代次数k要多大才使近似解x(k)的误差小于10-3？取x(0)=0.输入矩阵A，向量b，初值x0，eps计算B，g计算x1=B*x0+gwhilenorm(x1-x0,inf)>epsx0=x1k=k+1x1=B*x0+g输出k，x182/109解：在Matlab中B=[0.00000.1000-0.20000.0000;0.09090.00000.0909-0.2727;-0.20000.10000.00000.1000;0.0000-0.37500.12500.0000];g=[0.6000;2.2727;-1.1000;1.8750];norm(B,inf);x0=[0000]';x1=B*x0+g;k=1;whilenorm(x1-x0,inf)/norm(x1,inf)>0.001x0=x1;k=k+1x1=B*x0+gend对于先验估计，经计算需k=1383/109在简单迭代法中若将系数矩阵分裂为A=D+L+R，其中D为A的对角线部分，L为A的严格下三角部分，R为A的严格上三角部分。则可将方程组化为等价方程组：(D+L+R)x=bDx=(–L–R)x+b

x=D-1(–L–R)x+D-1b，从而得到迭代公式x(k+1)=D-1(–L–R)x(k)+D-1b，称为Jacobi迭代法，其迭代矩阵为：BJ=D-1(–L–R)。其分量式为84/109Jacobi迭代算法流程图：其中：BJ=D-1(–L–R)，g=D-1b输入矩阵A，向量b，初值x0，eps计算D，L，R计算B，g计算x1=B*x0+gwhilenorm(x1-x0,inf)>epsx0=x1k=k+1x1=B*x0+g输出k，x185/109例3-1用Jacobi迭代法解线性方程组取x(0)=0解：这里，，86/109Jacobi迭代法：87/109在Matlab中：A=[5-1-1-1;

-110-1-1;-1-15-1;-1-1-110];b=[-4;12;8;34];D0=diag(A);D=diag(D0);R=triu(A,1);L=tril(A,-1)BJ=inv(D)*(-L-R);g=inv(D)*bnorm(BJ,inf);x0=[0000]';x1=BJ*x0+g;k=1;whilenorm(x1-x0,inf)>0.001x0=x1;k=k+1x1=BJ*x0+gend88/109Jacobi迭代：x(k+1)=D-1(–L–R)x(k)+D-1b

x(k+1)=–D-1L

x(k)–D-1R

x(k)+D-1bSeidel迭代：x(k+1)=–D-1L

x(k+1)–D-1R

x(k)+D-1bDx(k+1)=–L

x(k+1)–R

x(k)+bDx(k+1)+L

x(k+1)=

–R

x(k)+bx(k+1)=

–(D+L)-1R

x(k)+(D+L)-1b即Seidel迭代的迭代矩阵：BS=-(D+L)-1R其分量式为89/109Seidel迭代的迭代矩阵：BS=-(D+L)-1R，g=(D+L)-1b算法：输入矩阵A，向量b，初值x0，eps计算D，L，R计算B，g计算x1=B*x0+gwhilenorm(x1-x0,inf)>epsx0=x1k=k+1x1=B*x0+g输出k，x190/109例3-2用Seidel迭代法解线性方程组取x(0)=0解：这里，，91/109Seidel迭代法：BS=-(D+L)-1R92/109在Matlab中：A=[5-1-1-1;

-110-1-1;-1-15-1;-1-1-110];b=[-4;12;8;34];D0=diag(A);D=diag(D0);R=triu(A,1);L=tril(A,-1);BS=-inv(D+L)*R;g=inv(D+L)*b;norm(BS,inf)x0=[0000]';x1=BS*x0+g;k=1;whilenorm(x1-x0,inf)>0.001x0=x1;k=k+1x1=BS*x0+gend93/1094.对角占优方程组定义：设A=(aij)nn，若，i=1，2，…，n，则称A按行严格对角占优

注：若A按行严格对角占优，则A可逆94/109定理：设Ax=b，若A按行严格对角占优，则对于任意的x(0)，J-迭代、S迭代均收敛。证明：(1)

，所以J-迭代收敛

95/109定理：设Ax=b，若A按行严格对角占优，则对于任意的x(0)，J-迭代、S迭代均收敛。证明：(2)BS=-(D+L)-1R——欲证(BS)<1因为：I–BS=I+(D+L)-1R=(D+L)-1[(D+L)+R]所以：0=|I–BS|=|(D+L)-1||

(D+L)+R||(D+L)+R|=0令：有|G|=0.若|

|1，则

，i=1，2，...，n

G严格对角占优

G可逆矛盾96/109注：若将方程组经初等变换化为同解方程组，使其系数矩阵严格对角占优，则可使迭代法收敛——p216例定理：若A正定，则S-迭代收敛若2D–A也正定，则J-迭代收敛，否则J-迭代不收敛97/109§2.5逐次超松弛迭代法和块迭代法1.逐次超松弛迭代法使用迭代法的困难是计算量难以估计，有些方程组的迭代格式虽然收敛，但收敛速度慢而使计算量变得很大。

逐次超松弛迭代法(SuccessiveOverRelaxationMethod,简写为SOR)可以看作带参数ω的高斯-塞德尔迭代法，是G-S方法的一种修正或加速.是求解大型稀疏矩阵方程组的有效方法之一.

这种方法将前一步的结果xi