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

1、第六章 线性方程组的直接法2/1092.0 引引 言言2.1 Gauss消去法消去法2.2 矩阵分解在解线性方程组中的作用矩阵分解在解线性方程组中的作用 补充补充 向量的范数与矩阵的范数向量的范数与矩阵的范数2.3 直接法的误差分析直接法的误差分析2.4 线性方程组迭代解法线性方程组迭代解法 3/1092.0 引引 言言 在自然科学和工程技术中很多问题的解决常常归结为在自然科学和工程技术中很多问题的解决常常归结为解线性代数方程组。解线性代数方程组。 例如电学中的网络问题，用最小二乘法求实验数据的例如电学中的网络问题，用最小二乘法求实验数据的曲线拟合问题，解非线性方程组问题，用差分法或者有限曲线

2、拟合问题，解非线性方程组问题，用差分法或者有限元方法解常微分方程、元方法解常微分方程、偏微分方程边值问题等都导致求解偏微分方程边值问题等都导致求解线性代数方程组。线性代数方程组。 而这些方程组的系数矩阵大致分为两种，一种是低阶而这些方程组的系数矩阵大致分为两种，一种是低阶稠稠密矩阵密矩阵(例如，阶数大约为例如，阶数大约为150)，另一种是大型稀疏矩，另一种是大型稀疏矩阵阵(即矩阵阶数高且零元素较多即矩阵阶数高且零元素较多)。 4/1092.0 引引 言言 设有线性方程组设有线性方程组Ax = = b，其中，其中 为非奇异阵，为非奇异阵， ， 关于线性方程组的数值解法一般有两类：直接法与迭关于线

3、性方程组的数值解法一般有两类：直接法与迭代法。代法。nnnnnnaaaaaaaaaA212222111211nxxxx21nbbbb215/1091. 直接法直接法 就是经过有限步算术运算，可求得方程组精确解的方就是经过有限步算术运算，可求得方程组精确解的方法法( (若计算过程中没有舍入误差若计算过程中没有舍入误差) )。 但实际计算中由于舍入误差的存在和影响，这种方法但实际计算中由于舍入误差的存在和影响，这种方法也只能求得线性方程组的近似解。也只能求得线性方程组的近似解。 本章将阐述这类算法中最基本的高斯消去法及其某些本章将阐述这类算法中最基本的高斯消去法及其某些变形。变形。 这类方法是解低

4、阶稠密矩阵方程组的有效方法，近十这类方法是解低阶稠密矩阵方程组的有效方法，近十几年来直接法在求解具有较大型稀疏矩阵方程组方面取得几年来直接法在求解具有较大型稀疏矩阵方程组方面取得了较大进展。了较大进展。6/1092. 迭代法迭代法 迭代法是用某种极限过程去逐步逼近线性方程组精确迭代法是用某种极限过程去逐步逼近线性方程组精确解的方法。解的方法。 迭代法具有需要计算机的存贮单元较少、程序设计简迭代法具有需要计算机的存贮单元较少、程序设计简单、原始系数矩阵在计算过程中始终不变等优点，但存在单、原始系数矩阵在计算过程中始终不变等优点，但存在收敛性及收敛速度问题。收敛性及收敛速度问题。 迭代法是解大型稀

5、疏矩阵方程组迭代法是解大型稀疏矩阵方程组( (尤其是由微分方程尤其是由微分方程离散后得到的大型方程组离散后得到的大型方程组) )的重要方法。的重要方法。7/1092.1 Gauss消去法消去法 高斯（高斯（GaussGauss）消去法是解线性方程组最常用的方法）消去法是解线性方程组最常用的方法之一。之一。 它的基本思想是通过逐步消元（行的初等变换），把它的基本思想是通过逐步消元（行的初等变换），把方程组化为系数矩阵为三角形矩阵的同解方程组，然后用方程组化为系数矩阵为三角形矩阵的同解方程组，然后用回代法解此三角形方程组（简单形式）得原方程组的解。回代法解此三角形方程组（简单形式）得原方程组的解。

6、 8/1092.1 Gauss消去法消去法 例：例：下面我们来讨论一般的解下面我们来讨论一般的解n阶方程组的高斯消去法。阶方程组的高斯消去法。21567003201111156140320111116122231111|bA9/1091. 消元消元将原方程组记为将原方程组记为A(1)(1)x = = b(1)(1)，其中其中A(1)(1)=(=(aij(1)(1) = () = (aij) )，b(1)(1) = = b(1) (1) 第一次消元。第一次消元。其中其中注：若注：若a1111(1)(1) = 0 = 0，则在第，则在第1 1列中至少有一个元素不为列中至少有一个元素不为0 0，可交

7、换行后再消元可交换行后再消元|00|)2()2()2()2(2)1(1)2()2(2)2(2)2(22)1(1)1(12)1(11)1()1(2)1(1)1()1(2)1(1)1(2)1(22)1(21)1(1)1(12)1(11)1()1(bAbbbaaaaaaabbbaaaaaaaaabAnnnnnnnnnnnnnniaabbaabbaanjaaaaaiiiiiijiijij,.,3 , 21,.,3 , 2)1(11)1(1)1(1)1(1)1(11)1(1)1()2()1(11)1(1)1()1(11)1(1)1()2(倍的减去倍行的减去第10/109(2) (2) 第第k次消元。次消

8、元。其中其中注：为减少计算量，令注：为减少计算量，令 ，则则|000|)1()1()1()1(1)()1(1)1()1(1,)1(, 1)1(1, 1)()(1,)()1(1)1(11)1(1)1(11)()()1(1)()()()()1(1)1(1)1(11)()(kkknkkkkknnkknknkkkkkknkkkkkknkkknkkknnknkkknkkknkkkbAbbbbaaaaaaaaaaabbbaaaaaaabAnkkiaabbaabbaaknkkjaaaaakkkkikkkkkkkkkikkikikkkkikkijkkkkikkijkij,.,2, 1,.,2, 1)()()(

9、)()()()()1()()()()()()()1(倍的减去倍行的减去第)()(kkkkikikaal nkkiblbbnkkjalaakkikkikikijikkijkij,.,2, 1,.,2, 1)()()1()()()1(11/109(3) (3) 当当k = = n 1 1时得时得完成第完成第n 1 1次消元后得到与原方程组等价的三角形方程次消元后得到与原方程组等价的三角形方程组组A( (n) )x = = b( (n) )注：当注：当det(A) 0时，显然有时，显然有aii(i) 0，(i = 1，n)，称，称为主元素。为主元素。)()2(2)1(1)()2(2)2(22)1(1

10、)1(12)1(11)()(|nnnnnnnnnbbbaaaaaabA12/1092. 回代回代求解三角形方程组求解三角形方程组A( (n) )x = = b( (n) )，得到求解公式，得到求解公式注：求解过程称为回代过程。注：求解过程称为回代过程。1,.,2, 1)(1)()()()(nniaxabxabxiiinijjiijiiinnnnnn13/1093. Gauss消去法的计算量消去法的计算量以乘除法的次数为主以乘除法的次数为主(1) 消元过程：消元过程：第第k步时步时(n k) + (n k) (n k + 1) = (n k) (n k + 2)共有共有nkkiblbbnkkja

11、laakkikkikikijikkijkij,.,2, 1,.,2, 1)()()1()()()1(nnnknknnk6523)2)(231114/1093. Gauss消去法的计算量消去法的计算量(2) 回代过程：回代过程：求求xi中，乘中，乘ni次，除次，除1 1次，共次，共ni+1+1次（次（i = 1，n1）共有共有1,.,2, 1)(1)()()()(nniaxabxabxiiinijjiijiiinnnnnn22) 1() 1(12111nnininnini15/1093. Gauss消去法的计算量消去法的计算量(3) (3) 总次数为总次数为注：当注：当n = 20 = 20时约

12、为时约为26702670次，比克莱姆法则次，比克莱姆法则9.79.7 10102020次大次大大减少。大减少。)(3333323nOnnnnMD16/1094. 说明说明(1) (1) 若消元过程中出现若消元过程中出现akk( (k) ) = 0 = 0，则无法继续，则无法继续(2) (2) 若若akk( (k) ) 0 0，但较小，则小主元做除数将产生大误，但较小，则小主元做除数将产生大误差差(3) (3) 每次消元都选择绝对值最大者作主元，称为高斯主元每次消元都选择绝对值最大者作主元，称为高斯主元消去法消去法(4) (4) 通常第通常第k步选择步选择第第k列主对角元以下元素绝对值最大者作主

13、元（该行列主对角元以下元素绝对值最大者作主元（该行与第与第k行对调），称为列主元消去法。行对调），称为列主元消去法。|max)(kiknikkac17/109例例1 1 用舍入三位有效数字求解线性方程组用舍入三位有效数字求解线性方程组（准确解是（准确解是x1 = 10.0，x2 = 1.00）解：解：1) 不选主元的不选主元的GaussGauss消去法计算结果：消去法计算结果： x1 = -10.0，x2 = 1.01，此解无效；此解无效； 2) 按列选主元的按列选主元的GaussGauss消去法计算结果：消去法计算结果： x1 = 10.0，x2 = 1.00. .0 .4710. 631.

14、 52 .599 .580300. 02121xxxx18/109(5) (5) 行标度化行标度化 当线性方程组的系数矩阵的元素大小相差很大时，则当线性方程组的系数矩阵的元素大小相差很大时，则可能引进因丢失有效数字而产生的误差，并且舍入误差的可能引进因丢失有效数字而产生的误差，并且舍入误差的影响严重，即使用影响严重，即使用GaussGauss主元消去法得到的解也不可靠主元消去法得到的解也不可靠 为避免这一问题，可将系数矩阵的行元素按比例增减为避免这一问题，可将系数矩阵的行元素按比例增减以使元素的变化减小以使元素的变化减小 如用每行元素的最大模除该行各元素，使它们的模都如用每行元素的最大模除该行

15、各元素，使它们的模都不大于不大于1 1，这叫行标度化，其目的是要找到真正的主元，这叫行标度化，其目的是要找到真正的主元19/109(5) (5) 行标度化行标度化 如用每行元素的最大模除该行各元素，使它们的模都如用每行元素的最大模除该行各元素，使它们的模都不大于不大于1 1，这叫行标度化，其目的是要找到真正的主元，这叫行标度化，其目的是要找到真正的主元 消元过程仍是对原方程组进行的，只不过在消元过程仍是对原方程组进行的，只不过在GaussGauss列列主元消去法的算法中，按列选主元主元消去法的算法中，按列选主元ck时，应修改为时，应修改为其中其中)()(|maxkikiknikksac.,.,

16、1,|,|max)()(nkkiaskijnjkki20/1095. 算法设计算法设计输入A，n，epsFor (k = 1 to n - 1)选主元A(I0，k) 确定行号p = A(I0，k)若| p | = epsText = 1，退出循环若I0 kT换行（I0 k）消元若| A(n，n) | | A(I0，k)|TI0 = ifor (i = k+1 to n)A(i，k) = A(i，k)/ A(k，k)for (j = k+1 to n+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

17、1 to 1 step -1)w = 0for (j = k+1 to n)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 = k to n)n = A(k，j), A(k，j) = A(I0，j), A(I0，j) = n21/1095. 算法设计算法设计(1) 高斯消去法高斯消去法a=5 -1 -1 -1 -4; -1 10 -1 -1 12; -1 -1 5 -1 8; -1 -1 -1 10 34;x=0;0;0;0n=4for k=1:n-1 for i=k+1:n a(i,:

18、)=a(i,:)-a(k,:)*a(i,k)/a(k,k) endendfor j=n:-1:1 x(j)=(a(j,n+1)-a(j,j+1:n)*x(j+1:n)/a(j,j)end22/109(2) (2) 列主元素消去法列主元素消去法a=-0.04 0.04 0.12 3;a=-0.04 0.04 0.12 3;0.56 -1.56 0.32 1;-0.24 1.24 -0.28 0;0.56 -1.56 0.32 1;-0.24 1.24 -0.28 0;x=0;0;0; n=3;x=0;0;0; n=3;for k=1:n-1for k=1:n-1 c,i=max(abs(a(k:

19、n,k); c,i=max(abs(a(k:n,k); q=i+k-1q=i+k-1 if q=k if q=k m=a(q,:);a(q,:)=a(k,:);a(k,:)=m m=a(q,:);a(q,:)=a(k,:);a(k,:)=m end end for i=k+1:n for i=k+1:n a(i,:)=a(i,:)-a(k,:) a(i,:)=a(i,:)-a(k,:)* *a(i,k)/a(k,k)a(i,k)/a(k,k) end endendendfor j=n:-1:1for j=n:-1:1 x(j)=(a(j,n+1)-a(j,j+1:n)x(j)=(a(j,n+1

20、)-a(j,j+1:n)* *x(j+1:n)/a(j,j)x(j+1:n)/a(j,j)endend23/1092.2 矩阵分解在解线性方程组中的作用矩阵分解在解线性方程组中的作用1. 原理原理若若A = = LR，则则Ax = = b LRx = = b 若其中若其中L、R为三角形矩阵，则方程组易解为三角形矩阵，则方程组易解定义：定义：(1) (1) 称称 为单位下三角阵为单位下三角阵(2) (2) 设设L为单位下三角阵，为单位下三角阵，R为上三角阵，若为上三角阵，若A = = LR，则，则称称A可三角（可三角（LR）分解，并称）分解，并称A = = LR为为A的三角分解（杜的三角分解（杜

21、利特尔分解）利特尔分解）yRxbLy1*1*1L24/109定理定理：(1) (1) A = ( = (aij) )n n有唯一有唯一LR分解分解 A的顺序主子的顺序主子式式 k 0 0，k = 1 = 1，2 2，.，n 1 1(2) (2) 若若A = ( = (aij) )n n可逆，则存在置换阵可逆，则存在置换阵P，使，使PA的的n个个顺顺序主子式全不等于序主子式全不等于0 0注：由注：由Ax = = b PAx = = Pb知，经适当行交换后，知，经适当行交换后，A总总可三角分解可三角分解25/1092. 直接三角分解法直接三角分解法设设A的各阶顺序主子式不为的各阶顺序主子式不为0，

22、且，且A = LR，即，即(1) (1) 主对角线（含）上边：当列标主对角线（含）上边：当列标m 行标行标k时，时，lkm = 0 = 0 j = = k，k+1+1，.，n；k = 1 = 1，2 2，.，n(2) (2) 主对角线下边：当行标主对角线下边：当行标m 列标列标k时，时，rkm = 0 = 0 i = = k+1+1，k+2+2，.，n；k = = 1，2，.，n1 1nnnnnnnnnnnnrrrrrrlllaaaaaaaaa22211211212121222211121111126/109设设A的各阶顺序主子式不为的各阶顺序主子式不为0，且，且A = LR，即，即(1) (

23、1) 主对角线（含）上边：当列标主对角线（含）上边：当列标m 行标行标k时，时，lkm = 0 = 0 j= =k,k+1, +1, ., ,n；k=1,2,=1,2,., ,n(2) (2) 主对角线下边：当行标主对角线下边：当行标m 列标列标k时，时，rkm = 0 = 0 i= =k+1,+1,k+2,+2,., ,n；k= =1,2,.,n1 1欲求欲求lik和和rkjnnnnnnnnnnnnrrrrrrlllaaaaaaaaa222112112121212222111211111kmmjkmnmmjkmkjrlrla11kmmkimnmmkimikrlrla1127/109(1) (

24、1) 主对角线（含）上边：当列标主对角线（含）上边：当列标m 行标行标k时，时，lkm = 0 = 0 j= =k,k+1, +1, ., ,n；k=1,2,=1,2,., ,n(2) (2) 主对角线下边：当行标主对角线下边：当行标m 列标列标k时，时，rkm = 0 = 0 i= =k+1,+1,k+2,+2,., ,n；k= =1,2,.,n1 1欲求欲求lik和和rkj设设k = 1 = 1，a1 1j = = l1111r1 1j r1 1j = = a1 1j j = 1 = 1，2 2，.，nai1 1 = = li1 1r1111 li1 1 = = ai1 1/ /r1111

25、 i = 2 = 2，3 3，.，nkmmjkmnmmjkmkjrlrla11kmmkimnmmkimikrlrla1128/109一般地，一般地，最后：最后：1,.3 , 2,.,2, 1,.,1,1111nknkkirrlalnkkjrlarkkkmmkimikikkmmjkmkjkj11nmmnnmnnnnrlarr11r12.r1n第1步l21r22.r2n第2步l31l32.ln1ln2rnn第n步29/1093. 求解方程组求解方程组下三角下三角 Ly = = b：上三角上三角 Rx = = y：紧凑格式：紧凑格式：niylbybyijjijii,.,3 , 211111,.,2,

26、 1/ /1nnirxryxryxiinijjijiinnnnnnnnnnnnnnnnnnnnnnyyyyrlllrrllrrrlrrrrbbbbaaaaaaaaaaaaaaaa32132133332312232221113121132132133332312232221113121130/109a=2 10 0 -3; -3 -4 -12 13; 1 2 3 -4; 4 14 9 -13;a=2 10 0 -3; -3 -4 -12 13; 1 2 3 -4; 4 14 9 -13;b=10;5;-2;7; y=0 0 0 0b=10;5;-2;7; y=0 0 0 0; x=0 0 0 0

27、; x=0 0 0 0n=4;n=4;for r=1:nfor r=1:n a(r,r:n) = a(r,r:n)-a(r,1:r-1) a(r,r:n) = a(r,r:n)-a(r,1:r-1)* *a(1:r-1,r:n)a(1:r-1,r:n) a(r+1:n,r)=(a(r+1:n,r)-a(r+1:n,1:r-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)a(1:r-1,r)/a(r,r)endendR = triu(a,0)R = triu(a,0)L =eye(n)+ tril(a,-1)L =eye(n

28、)+ tril(a,-1)for r=1:nfor r=1:n y(r) = b(r)-L(r,1:r-1) y(r) = b(r)-L(r,1:r-1)* *y(1:r-1)y(1:r-1)endendfor j=n:-1:1for j=n:-1:1 x(j)=(y(j)-R(j,j+1:n) x(j)=(y(j)-R(j,j+1:n)* *x(j+1:n)/R(j,j)x(j+1:n)/R(j,j)endend31/109紧凑格式紧凑格式a=2 10 0 -3 10;a=2 10 0 -3 10; -3 -4 -12 13 5; -3 -4 -12 13 5; 1 2 3 -4 -2; 1

29、 2 3 -4 -2; 4 14 9 -13 7; 4 14 9 -13 7;x=0 0 0 0 x=0 0 0 0n=4;n=4;for r=1:nfor r=1:n a(r,r:n+1) = a(r,r:n+1)-a(r,1:r-1) a(r,r:n+1) = a(r,r:n+1)-a(r,1:r-1)* *a(1:r-1,r:n+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(r+1:n,r)=(a(r+1:n,r)-a(r+1:n,1:r-1)* *a(1:r-1,r)/a(r,r)a(1:r-1,r)/a(r,r)en

30、dendR = triu(a,0)R = triu(a,0)for j=n:-1:1for j=n:-1:1 x(j)=(a(j,n+1)-R(j,j+1:n) x(j)=(a(j,n+1)-R(j,j+1:n)* *x(j+1:n)/R(j,j)x(j+1:n)/R(j,j)endend32/1094. 选主元的三角分解法选主元的三角分解法5. 平方根法平方根法(1) (1) 定理定理：若若A正定正定，则则A可唯一分解为可唯一分解为A = = LDLT T。其中。其中L为单位下三角，为单位下三角，D为元素全为正数的对角阵。为元素全为正数的对角阵。(2) (2) 设设A = = LDLT T则

31、有：则有：1111112121212121TnnnnnllldddlllLDLAnknkkiddllaldladkkmmkmimikikkmmkmkkk,.,2 , 1,.,2, 11111233/109(3) (3) 为了避免重复计算，引进中间量为了避免重复计算，引进中间量tik = = likdk，将上式改，将上式改写为：写为：nknkkidtlltatltadkikikkmkmimikikkmkmkmkkk,.,2 , 1,.,2, 1/1111nknkkiddllaldladkkmmkmimikikkmmkmkkk,.,2 , 1,.,2, 11111234/109(4) (4) Ax

32、 = = b LDLT Tx = = b Ly = = b，Dz = = y，LT Tx = = z求解公式：求解公式：上述方法称为上述方法称为“改进的平方根法改进的平方根法”注：不能选主元作行交换注：不能选主元作行交换破坏对称性，必是数值稳定破坏对称性，必是数值稳定的的) 1,.,1,(),.,2 , 1(/),.,2 , 1(111nnixlzxnidyzniylbynimmmiiiiiiimmimii35/109a=5 -4 1; -4 6 -4; 1 -4 6; b=2;-1;-1a=5 -4 1; -4 6 -4; 1 -4 6; b=2;-1;-1d=0 0 0; t=0 0 0;

33、0 0 0;0 0 0; L=0 0 0;0 0 0;0 0 0;d=0 0 0; t=0 0 0;0 0 0;0 0 0; L=0 0 0;0 0 0;0 0 0;n=3n=3for k=1:nfor k=1:n d(k)=a(k,k)- t(k,1:k-1)d(k)=a(k,k)- t(k,1:k-1)* *L(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)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:k-1) L(k+1:n,k)=t(k+1:n,k)/d(k)

34、 L(k+1:n,k)=t(k+1:n,k)/d(k)endendx=0 0 0;y=0 0 0;z=0 0 0;x=0 0 0;y=0 0 0;z=0 0 0;for i=1:nfor i=1:n y(i)=b(i)-L(i,1:i-1) y(i)=b(i)-L(i,1:i-1)* *y(1:i-1)y(1:i-1)endendfor i=1:nfor i=1:n z(i)=y(i)/d(i) z(i)=y(i)/d(i)endendfor i=n:-1:1for i=n:-1:1 x(i)=z(i)-L(i+1:n,i) x(i)=z(i)-L(i+1:n,i)* *x(i+1:n)x(i

35、+1:n)endend36/1096. 追赶法追赶法 在实际问题中，经常遇到以下形式的方程组在实际问题中，经常遇到以下形式的方程组这种方程组的系数矩阵称为三对角矩阵。以下针对该方程这种方程组的系数矩阵称为三对角矩阵。以下针对该方程组的特点提供一种简便有效的算法组的特点提供一种简便有效的算法 追赶法。追赶法。nnnnnnnnnddddxxxxbacbacbacbdAx1211211112221137/109设直接三角分解为：设直接三角分解为：A = = LR其中其中容易看出容易看出pi = = ci（1 1，2 2，n 1 1），计算），计算L和和R的其余的其余元素的公式为元素的公式为,1112

36、nllLnnrprprR1211),.,3 , 2(1111niclbrralbriiiiiii,12211nnnbacbacbA38/109容易看出容易看出pi = = ci（1 1，2 2，n 1 1），计算），计算L和和R的其余的其余元素的公式为元素的公式为计算过程：计算过程：r1 1 l2 2 r2 2 . . ln rn),.,3 , 2(1111niclbrralbriiiiiii39/109有了有了A的的LR分解后，线性方程组分解后，线性方程组Ax = = d等价于两个简单的等价于两个简单的方程组：方程组：Ly = = d，Rx = = y于是，计算于是，计算y的公式是的公式是y

37、1 1 = = d1 1，yi = = di liyi1 1 （i = 2，3，n）计算过程：计算过程：r1 1y1 1l2 2r2 2y2 2.lnrnyn （追）（追）计算计算x的公式是的公式是xn= =yn/ /rn，xi=(=(yicixi+1+1)/)/ri （i=n1，n2，1）( (赶）赶）nnnnnnnyyyxxxrcrcrdddyyyll2121121121212,11140/109定理：若定理：若A为对角占优（为对角占优（diagonally dominantdiagonally dominant）的三对）的三对角阵，且满足角阵，且满足|b1| |c1| 0，|bi| |a

38、i| + |ci| ，aici 0，i =2，3，n-1|bn| |an| 0.则追赶法可解以则追赶法可解以A为系数矩阵的方程组为系数矩阵的方程组注：注：1. 1. 如果如果A是严格对角占优阵，则不要求三对角线上的所是严格对角占优阵，则不要求三对角线上的所有元素非零。有元素非零。2. 2. 根据不等式根据不等式可知：分解过程中，矩阵元素不会过分增大，算法保证稳可知：分解过程中，矩阵元素不会过分增大，算法保证稳定。定。3. 3. 运算量为运算量为 O(6O(6n) )。|,1|/|iiiiiiiabrabrc41/1097. QR算法算法42/109补充补充 向量的范数与矩阵的范数向量的范数与矩

39、阵的范数 在线性方程组的数值解法中，经常需要分析解向量的在线性方程组的数值解法中，经常需要分析解向量的误差，需要比较误差向量的误差，需要比较误差向量的“大小大小”或或“长度长度”。那么怎。那么怎样定义向量的长度呢？样定义向量的长度呢？ 我们在初等数学里知道，定义向量的长度，实际上就我们在初等数学里知道，定义向量的长度，实际上就是对每一个向量按一定的法则规定一个非负实数与之对应，是对每一个向量按一定的法则规定一个非负实数与之对应，这一思想推广到维线性空间里，就是向量的范数或模。这一思想推广到维线性空间里，就是向量的范数或模。 用用Rn表示表示n维实向量空间，用维实向量空间，用Cn表示表示n维复向

40、量空间，维复向量空间，首先将向量长度概念推广到首先将向量长度概念推广到Rn( (或或Cn) )中。中。43/1091. 向量的范数向量的范数向量的范数可以看作是描述向量向量的范数可以看作是描述向量“大小大小”的一种度量的一种度量. .范数的最简单的例子，是绝对值函数：范数的最简单的例子，是绝对值函数：有三个熟知的性质：有三个熟知的性质：(1) x 0 x 0 x = 0当且仅当当且仅当x = 0(2) ax = a x a为常数为常数(3) x + y x + y 2xx 44/1091. 向量的范数向量的范数范数的另一个简单例子是三维欧氏空间的长度范数的另一个简单例子是三维欧氏空间的长度设设

41、x = ( = (x1 1, , x2 2, , x3 3) )，则，则x的欧氏范数定义为：的欧氏范数定义为：欧氏范数也满足三个条件：欧氏范数也满足三个条件： x，y R R3 3，a为常数为常数(1) (1) x 0 0，且且 x = 0 = 0 x = 0= 0(2) (2) ax = = a x (3) (3) x + + y x + + y 前两个条件显然，第三个条件在几何上解释为三角形一边前两个条件显然，第三个条件在几何上解释为三角形一边的长度不大于其它两边长度之和。因此，称为三角不等式。的长度不大于其它两边长度之和。因此，称为三角不等式。232221|xxxx45/109向量范数的

42、一般概念：向量范数的一般概念：定义定义1 1：设设V是数域是数域F上的向量空间，对上的向量空间，对V中任一向量中任一向量，都，都有唯一实数有唯一实数 与之对应，满足如下三个条件：与之对应，满足如下三个条件：1) 1) 正定性：正定性： 0 0，且，且 = 0 = 0 = 0= 02) 2) 齐次性：齐次性： k = | = |k| | ，这里，这里k F3) 3) 三角不等式：三角不等式： + + + + 则称则称 为为的范数。定义了范数的向量空间称为赋范向量的范数。定义了范数的向量空间称为赋范向量空间空间. .简单性质：简单性质：(1) (1) x 0 单位向量单位向量(2) |x| = |

43、 x |(3) | |x| |y| | | x y | 当x y时，|x| |y|1|xx46/109C Cn上的常见范数有：上的常见范数有：1) 1-1) 1-范数范数 2) 2-2) 2-范数范数 称为欧氏范数称为欧氏范数3) 3) - -范数范数不难验证，上述三种范数都满足定义的条件。不难验证，上述三种范数都满足定义的条件。注：上述形式的统一：注：上述形式的统一： 1 1 p niixx11|niixx122|max|1inixxpnipipxx/11)|(|47/109定理定理5 5：定义在定义在C Cn上的向量范数上的向量范数|x|是变量是变量x分量的一致连分量的一致连续函数。（续函

44、数。（f( (x) = ) = |x|）定理定理6 6：在在C Cn上定义的任何两个范数都是等价的。上定义的任何两个范数都是等价的。即存在正数即存在正数k1 1与与k2 2( (k1 1 k2 2 0)0)，对一切，对一切x Cn，不等式，不等式k1 1| x |b | x |a k2 2| x |b 成立。成立。对常用范数，容易验证下列不等式：对常用范数，容易验证下列不等式：111xxxnxnxx1xnxx248/109例例1 1 设设A为正定矩阵，为正定矩阵， x Cn，令，令 ，则则| x |A为向量范数，称为加权范数为向量范数，称为加权范数证明：证明：(1) (1) 显然显然 x A

45、0 0，且，且 x A = 0 = 0 x = 0= 0(2) (2) (3) (3) 因为因为A正定，所以正定，所以 可逆可逆P，使，使A = = PT TP 所以所以|x+ +y|A= =|P( (x+ +y) )|2 2= =|Px+ +Py|2 2 |Px|+ +|Py|2 2 = = |x|A+ +| y|A综上，综上，| x+ +y |A为范数。为范数。注：注：当当0 p 1时，时， 不是范数不是范数AxxxAT|AAxkAxxkAxxkkxAkxkx|)()(|TT2T2TTTT|)()(|PxPxPxPxPxAxxxApnipipxx/11)|(|49/109有了范数的概念，就

46、可以讨论向量序列的收敛性问题。有了范数的概念，就可以讨论向量序列的收敛性问题。定义定义2：设给定设给定Cn中的向量序列中的向量序列 xk ，即，即x0 0，x1 1，xk，其中其中若对任何若对任何i ( (i = 1, 2, = 1, 2, , n) )都有都有则向量则向量 称为向量序列称为向量序列 xk 的极限，或者说的极限，或者说向量序列向量序列 xk 依坐标收敛于向量依坐标收敛于向量x* *，记为，记为T)()(2)(1.,knkkkxxxx *)(limikikxxT*1*).,(nxxx *limxxkk50/109定义定义2：设给定设给定Cn中的向量序列中的向量序列 xk ，即，即

47、x0 0，x1 1，xk，其中其中若对任何若对任何i ( (i = 1, 2, = 1, 2, , n) )都有都有则向量则向量 称为向量序列称为向量序列 xk 的极限，或者说的极限，或者说向量序列向量序列 xk 依坐标收敛于向量依坐标收敛于向量x* *，记为，记为定理定理7 7：向量序列向量序列 xk 依坐标收敛于依坐标收敛于x* *的充要条件是的充要条件是如果一个向量序列如果一个向量序列 xk 与向量与向量x* *满足上式，就说向量序列满足上式，就说向量序列 xk 依范数收敛于依范数收敛于x* *，于是便得：向量序列依范数收敛与，于是便得：向量序列依范数收敛与依坐标收敛是等价的。依坐标收敛

48、是等价的。T)()(2)(1.,knkkkxxxx *)(limikikxxT*1*).,(nxxx *limxxkk0lim*xxkk51/1092. 矩阵的范数矩阵的范数下面我们给出矩阵范数的一般定义。为简单起见，仅讨论下面我们给出矩阵范数的一般定义。为简单起见，仅讨论实数域的情形。实数域的情形。定义定义3 3 若矩阵若矩阵A Rnn的某个非负的实值函数的某个非负的实值函数N( (A)= )= |A|，满足条件满足条件1) 1) 正定性正定性： A 0 0，且且 A = 0 = 0 A = 0= 02) 2) 齐次性齐次性： kA = = |k| |A| k R3) 3) 三角不等式三角不

49、等式：| A + + B | | A | + + | B |4) 4) 相容性相容性：|AB| | A | B |则称则称N( (A) )是是Rnn上的一个矩阵范数上的一个矩阵范数( (或模或模) )。52/1092. 矩阵的范数矩阵的范数定义定义3 3 若矩阵若矩阵A Rnn的某个非负的实值函数的某个非负的实值函数N( (A)= )= |A|，满足条件满足条件1) 1) 正定性正定性： A 0 0，且且 A = 0 = 0 A = 0= 02) 2) 齐次性齐次性： kA = = |k| |A| k R3) 3) 三角不等式三角不等式：| A + + B | | A | + + | B |4

50、) 4) 相容性相容性：|AB| | A | B |则称则称N( (A) )是是Rnn上的一个矩阵范数上的一个矩阵范数( (或模或模) )。如由如由Rnn上上2-范数可以得到范数可以得到Rnn中矩阵的一种范数中矩阵的一种范数称为称为A的的FrobeniusFrobenius范数。范数。|A|F显然显然满足正定性、齐次性满足正定性、齐次性及三角不等式。及三角不等式。2/11,2)(|)(njiijFaAAF53/109由于在大多数与估计有关的问题中，矩阵和向量会同时参由于在大多数与估计有关的问题中，矩阵和向量会同时参与讨论，所以希望引进一种矩阵的范数，它是和向量范数与讨论，所以希望引进一种矩阵的

51、范数，它是和向量范数相联系而且和向量范数相容的，即相联系而且和向量范数相容的，即5) 5) |Ax| |A| |x|对任何向量对任何向量x Rn及矩阵及矩阵A Rnn都成立。都成立。注：满足相容性条件的范数很多，其中最重要的是算子范注：满足相容性条件的范数很多，其中最重要的是算子范数数54/109再引进一种矩阵的范数。再引进一种矩阵的范数。定义定义4 4 设设x Rn及矩阵及矩阵A Rnn，|x|v为一种向量范数为一种向量范数(如如v = 1，2或或)，相应地定义一个矩阵，相应地定义一个矩阵的非负函数的非负函数可验证可验证|A|v满足定义满足定义3，所以，所以|A|v是是Rnn上矩阵的一个范上

52、矩阵的一个范数，数，|A|v还满足相容性条件还满足相容性条件5)，称为，称为A的算子范数，或诱的算子范数，或诱导范数。导范数。注注：(1) 可以证明：可以证明：(2) |Ax|是是x的连续函数，从而存在的连续函数，从而存在x*使使|x*| = 1且且|Ax*| = |A|(3) 若若| . |为算子范数，为算子范数，I为单位矩阵，则为单位矩阵，则|I| = 1vvxvxAxA|max|0vxvAxA|max|1|55/109定理定理8 8：设设n阶方阵阶方阵A = ( = (aij) )n n，则则与与|x|1相容相容的矩阵范数是的矩阵范数是与与|x|2相容相容的矩阵范数是的矩阵范数是其中其中

53、 1 1为矩阵为矩阵AT TA的最大特征值。的最大特征值。与与|x| 相容的矩阵范数是相容的矩阵范数是上述范数分别称为上述范数分别称为1-1-范数，范数，2-2-范数和范数和 -范数范数又称为列模、谱模和行模又称为列模、谱模和行模niijjaA11|max|12|AnjijiaA1|max|56/109注：注：A的范数的范数|A|与与A的特征值间的关系的特征值间的关系设设 为矩阵为矩阵A的任一特征值，的任一特征值，e为相应的特征向量，则为相应的特征向量，则 e= Ae，因为因为| | |e| = |Ae| |A| |e|，所以所以| | |A|从而得从而得定理定理9 9：矩阵矩阵A的任一特征值

54、的绝对值不超过的任一特征值的绝对值不超过A的范数的范数|A|。定义定义6 6：矩阵矩阵A的诸特征值的最大绝对值称为的诸特征值的最大绝对值称为A的谱半径，的谱半径，记为：记为：由定理由定理9 9可知：可知：iniA1max)(AA )(57/109例例(p138)(p138)：设：设A = (aij)n n，|A|为其算子范数，则为其算子范数，则|A| 1 I A可逆，且|11|)(|1AAI58/1092.3 直接法的误差分析直接法的误差分析1. 误差向量与残差向量误差向量与残差向量设设x为线性方程组为线性方程组Ax = = b的准确解，用数值算法得到的近的准确解，用数值算法得到的近似解为似解

55、为x* *，则称向量，则称向量e = = x* * - - x为误差向量。为误差向量。显然，显然，| e |可以表示近似解可以表示近似解x* *的准确程度：的准确程度：| e |越小，近越小，近似解似解x* *越准确。但由于越准确。但由于x不知道，所以无法计算不知道，所以无法计算e。变通的办法是考虑残差向量：变通的办法是考虑残差向量：r = = bAx* *的范数的大小。的范数的大小。但是，在某些情况下，尽管但是，在某些情况下，尽管| r |很小，很小，| e |也可能很大。也可能很大。59/109例例 线性方程组线性方程组的准确解是的准确解是(1(1，1)1)T T，若用某种方法得到计算解，

56、若用某种方法得到计算解(2.000(2.000，0.500)0.500)，则残差向量，则残差向量其其 - -范数范数| r | = 0.0015，而误差向量而误差向量其其 - -范数范数| e | = 1，它是残差向量的它是残差向量的667667倍。倍。500. 1000. 3001. 14990. 0000. 2000. 121xx0015. 00500. 1000. 3500. 0000. 2001. 14990. 0000. 2000. 1r500. 0000. 1e60/109能否建立残差向量与误差向量之间的关系？能否建立残差向量与误差向量之间的关系？由于由于r = b Ax* = A

57、x Ax* = A(x x*) = Ae，故有故有e = A1r，|e| |A1|r|另一方面另一方面，由由b = Ax得得|b| |A|x|，因此因此其中，其中， 是近似解是近似解x* *的相对误差，而的相对误差，而 表示相对残差表示相对残差上式表明近似解上式表明近似解x* *的相对误差大约是相对残差的的相对误差大约是相对残差的|A|A1|倍倍。|111brAAxArAAxrAxe|xe|br61/1092. 条件数与病态方程组条件数与病态方程组定义定义5 5：设：设A为为n阶非奇矩阵，称数阶非奇矩阵，称数|A|A1|为为矩阵矩阵A的条件的条件数，记为数，记为cond(A)。 显然，当显然，

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

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

60、也非易事，所以一般地说，遇到下列几种情况之一就应考也非易事，所以一般地说，遇到下列几种情况之一就应考虑线性方程组可能是病态的：虑线性方程组可能是病态的： (1) 系数矩阵的某两行系数矩阵的某两行( (列列) )几乎近似相关；几乎近似相关； (2) 系数矩阵的行列式的值很小；系数矩阵的行列式的值很小； (3) 用主元消去法解线性方程组时出现小主元；用主元消去法解线性方程组时出现小主元； (4) 近似解近似解x*已使残差向量已使残差向量r = b - Ax*的范数很小，但的范数很小，但该近似解仍不符合问题的要求该近似解仍不符合问题的要求64/1094. 解的精度估计与改善解的精度估计与改善 我们知