第三章 MATLAB线性方程组.ppt

1、第三章 线性代数方程组及矩阵特征值,3.0 预备知识 3.1 直接法 3.2 迭代法 3.3 不可解问题 3.4 病态问题, 3.0预备知识,一、对角阵与三角阵 1、对角阵： diag(A) 提取mn的矩阵A 的主对角线上元素，生 成一个具有min(m,n)个元素的列向量 diag(A,k) 提取第k条对角线的元素 diag(V) 设V为具有m个元素的向量，将产生一个以向量V的元素为主对角线元素的mm对角矩阵diag(V,k) 产生一个以向量V的元素为第k条对角线的元素的nn(n=m+k)对角阵,2、 矩阵的三角阵 下三角矩阵 tril(A) 提取矩阵A的下三角阵 tril(A,k) 提取矩阵

2、A的第k条对角线以下的元素 上三角矩阵 triu(A) 提取矩阵A的上三角矩阵 triu(A,k) 提取矩阵A的第k条对角线以下的元素,二、 用于专门学科中的矩阵 (1) 魔方矩阵 魔方矩阵是每行、每列及两条对角线上的元素和都相等的矩阵。对于n阶魔方阵，其元素由1,2,3,n2共n2 个整数组成. magic(n) 生成一个n阶魔方阵,各行各列及两条 对角线和为(1+2+3+.+n2 )/n,例 magic(5) ans = 17 24 1 8 15 23 5 7 14 16 4 6 13 20 22 10 12 19 21 3 11 18 25 2 9 例：将101125等25个数填入一个5

3、行5列的表格中，使其每行每列及对角线的和均为565。命令如下： M=100+magic(5),(2) 范得蒙矩阵 范得蒙(Vandermonde)矩阵最后一列全为1，倒数第二列为一个指定的向量，其他各列是其后列与倒数第二列的点乘积。 vander(V) 生成以向量V为基础向量的范得蒙矩阵。 例 A=vander(1;2;3;5) A = 1 1 1 1 8 4 2 1 27 9 3 1 125 25 5 1,(3) 希尔伯特矩阵（Hilbert） Hilbert矩阵的每个元素hij=1/(i+j-1) hilb(n) 生成n阶希尔伯特矩阵 invhilb(n) 专求n阶的希尔伯特矩阵的逆矩阵

4、注意1：高阶Hilbert矩阵一般为病态矩阵，所以直接求逆可能出现错误结论，故应该采用invhilb(n) 注意2：由于Hilbert矩阵本身接近奇异，所以建议处理该矩阵时建议尽量采用符号运算工具箱，若采用数值解时应该验证结果的正确性,(4) 托普利兹矩阵 (toeplitz) toeplitz矩阵除第一行第一列外，其他每个元素都与左上角的元素相同。 toeplitz(x,y) 生成一个以x为第一列，y为第一行的托普利兹矩阵。这里x, y均为向量，两者不必等长。toeplitz(x)用向量x生成一个对称的托普利兹矩阵。 例 T=toeplitz(1:5,1:7) T = 1 2 3 4 5 6

5、 7 2 1 2 3 4 5 6 3 2 1 2 3 4 5 4 3 2 1 2 3 4 5 4 3 2 1 2 3,(5) 帕斯卡矩阵 二次项(x+y)n展开后的系数随n的增大组成一个三角形表，称为杨辉三角形。由杨辉三角形表组成的矩阵称为帕斯卡(Pascal)矩阵。 pascal(n) 生成一个n阶帕斯卡矩阵。,pascal(6) ans = 1 1 1 1 1 1 1 2 3 4 5 6 1 3 6 10 15 21 1 4 10 20 35 56 1 5 15 35 70 126 1 6 21 56 126 252 例 求(x+y)5的展开式。 pascal(6)次对角线上的元素1，5，

6、10，10，5，1即为展开式的系数。,三、向量和矩阵的范数 norm(V)或norm(V,2) 求向量V（或矩阵V )的2范数 norm(V,1) 求向量V（或矩阵V）的1范数 norm(V,inf) 求向量V（或矩阵V)的范数 例 已知V，求V的3种范数。 命令如下： V=-1,1/2,1; v1=norm(V,1) %求V的1范数 v2=norm(V) %求V的2范数 vinf=norm(V,inf) %求V的范数,常用的向量范数:,范数意义下的单位向量: X=x1, x2T,-1,常用的矩阵范数:,例 求矩阵A的三种范数 命令如下： A=17,0,1,0,15;23,5,7,14,16;

7、4,0,13,0,22;10,12,19,21,3;11,18,25,2,19; a1=norm(A,1) %求A的1范数 a2=norm(A) %求A的2范数 ainf=norm(A,inf) %求A的范数,四、 矩阵的逆与伪逆 1、 矩阵的逆（后面研究完Gauss消去法后还将给出求逆的方法） 求方阵A的逆可用 inv(A) 注意：该函数也适用于符号变量构成的矩阵的求逆 例 用求逆矩阵的方法解线性方程组。 命令如下： A=1,2,3;1,4,9;1,8,27; b=5,2,6; x=inv(A)*b 一般情况下，用左除x=Ab比求矩阵的逆的方法更有效(因为A奇异或接近奇异时，用inv(A)可

8、能出错),例：Hilbert矩阵（非常有名的病态矩阵）：,验证从55到1414的Hilbert矩阵病态特征,clear format long; for n=5:14 H=hilb(n);norm1=norm(H*inv(H)-eye(size(H); H1=invhilb(n);norm2=norm(H*H1-eye(size(H); fprintf( n=%3.0f norm1=%e norm2=%en, n,norm1,norm2) end,n= 5 norm1=1.409442e-011 norm2=3.637979e-012 n= 6 norm1=2.534893e-009 norm

9、2=1.455203e-011 n= 7 norm1=9.810538e-008 norm2=5.208793e-010 n= 8 norm1=3.123671e-006 norm2=4.822187e-008 n= 9 norm1=8.116595e-005 norm2=2.479206e-007 n= 10 norm1=2.645008e-003 norm2=1.612897e-005 n= 11 norm1=7.200720e-002 norm2=1.305122e-003 Warning: Matrix is close to singular or badly scaled. Res

10、ults may be inaccurate. RCOND = 2.632766e-017. n= 12 norm1=1.176913e+001 norm2=5.576679e-002 Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.339949e-018. n= 13 norm1=5.323696e+001 norm2=1.137063e+001 Warning: Matrix is close to singular or badly scaled. Res

11、ults may be inaccurate. RCOND = 1.708191e-019. n= 14 norm1=1.224232e+004 norm2=2.836298e+002,说明1:对于Hilbert求逆时，不建议用inv,可用 invhilb直接产生逆矩阵 说明2：符号工具箱中也对符号矩阵定义了inv( )函数,即使对更高阶的奇异矩阵也可以精确的求解出逆矩阵,例：H=sym(hilb(30); norm(double(H*inv(H)-eye(size(H) 结果为 ans = 0,说明3：对于奇异矩阵用数值方法的inv( )函数，会给出错误的结果，但符号工具箱中的inv( )会

12、给出正确的结论,例 A=16,2,3,13;5,11,10,8; 9,7,6,12;4,14,15,1; det(A) B=inv(A),结果： ans = 0 Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.306145e-017. B = 1.0e+014 * 0.93824992236885 2.81474976710656 -2.81474976710656 -0.93824992236885 2.81474976710656 8.44424930131

13、968 -8.44424930131968 -2.81474976710656 -2.81474976710656 -8.44424930131968 8.44424930131968 2.81474976710656 -0.93824992236885 -2.81474976710656 2.81474976710656 0.93824992236885,但用符号工具箱的inv:,A=16,2,3,13; 5,11,10,8; 9,7,6,12; 4,14,15,1; A=sym(A) inv(A),结果： ? Error using = sym.inv Error, (in inverse

14、) singular matrix,2、 矩阵的伪逆 pinv(A) 若A不是一个方阵，或A是一个非满秩的方阵时，矩阵A没有逆矩阵，但可以找到一个与A的转置矩阵A同型的矩阵B，使得： ABA=A, BAB=B 此时称矩阵B为矩阵A的伪逆。 例 求A的伪逆，并将结果送B A=3,1,1,1;1,3,1,1;1,1,3,1; B=pinv(A) 例 求矩阵A的伪逆 A=0,0,0;0,1,0;0,0,1; pinv(A),五、 求方阵A的行列式: det(A) 例 用克莱姆(Cramer)方法求解线性方程组(不建议使用） 程序如下： D=2,2,-1,1;4,3,-1,2;8,5,-3,4;3,3

15、,-2,2; %定义系数矩阵 b=4;6;12;6; %定义常数项向量 D1=b,D(:,2:4); %用方程组的右端向量置换D的第1列 D2=D(:,1),b,D(:,3:4); %用方程组的右端向量置换D的第2列 D3=D(:,1:2),b,D(:,4:4);%用方程组的右端向量置换D的第3列 D4=D(:,1:3),b;%用方程组的右端向量置换D的第4列,DD=det(D); x1=det(D1)/DD; x2=det(D2)/DD; x3=det(D3)/DD; x4=det(D4)/DD; x1,x2,x3,x4,六、 求矩阵A的秩 rank(A) 七、 求矩阵A的迹 trace(A

16、) 例：D=2,2,-1,1;4,3,-1,2;8,5,-3,4;3,3,-2,2; r=rank(D) t=trace(A) 结果：r = 4 t= 6,九、求矩阵特征多项式、特征值、特征向量,设A是一个 nn 矩阵，,为A的特征多项式。,特征多项式的根称为矩阵A的特征值。,c=poly(A) 求矩阵A的特征多项式的系数 roots(c) 求多项式c的根,八、 求矩阵A的共轭矩阵 conj(A),eig(A) 求矩阵A的特征值 常用的调用格式有 ： E=eig(A) 求矩阵A的全部特征值，构成向量E。 V,D=eig(A) 求矩阵A的全部特征值，构成对角阵D，并求A的特征向量构成V的列向量。

17、 V,D=eig(A,nobalance) 与第2种格式类似，但第2种格式中先对A作相似变换后求矩阵A的特征值和特征向量，而格式3直接求矩阵A的特征值和特征向量，即A中某项非常小，这样求出的特征值及特征向量更精确。 V,D=eig(A,B) 计算广义特征值和特征向量，使 AV=BVD。,例：设矩阵,A=3 4 -2;3 -1 1;2 0 5; E=eig(A) V,D=eig(A) V,D=eig(A,nobalance),现求解线性方程组,情形1：m=n（恰定方程组） 在MATLAB中的求解命令有：,情形2：mn（超定方程），多用于曲线拟合。,解线性方程组的一般函数文件如下： functio

18、n x,y=line_solution(A,b) m,n=size(A);y=; if norm(b)0 % b非零，非齐次方程组 if rank(A)=rank(A,b) %方程组相容 （有解） if rank(A)=n %有唯一解 x=Ab; else %方程组有无穷多个解,基础解系 disp(原方程组有有无穷个解，其齐次方程组的基础 解系为y，特解为x); y=null(A,r); x=Ab; end,else %方程组不相容（无解） ，给出最小二乘解 disp(方程组的最小二乘法解是：); x=Ab; end else %齐次方程组 if rank(A)=n %列满秩 x=zero(n

19、,1) %0解 else %非0解，无穷多个解 disp(方程组有无穷个解，基础解系为x); x=null(A,r); end end return,如在MATLAB命令窗口，输入命令 A=2,2,-1,1;4,3,-1,2;8,5,-3,4;3,3,-2,2;b=4,6,12,6; x,y=line_solution(A,b) 及： A=2,7,3,1;3,5,2,2;9,4,1,7;b=6,4,2; x,y=line_solution(A,b) 分别显示其求解结果。,求解线性方程组（m=n),用克莱姆法则，理论上可行，但实际运算时工作量大，耗时。下面研究,、 一、Gauss简单消元法（m=

20、n),设 有,用线性代数中的克莱姆法则求解时，工作量非常大， 工作量小的方法是 Gauss消元法。,3.1 解线性方程组的直接法,消 元：,以此类推，最后方程组化为：,回 代：,失效，故应选主元,二、列主元素消去法-计算结果可靠,到此原方程组化为,到此原方程组化为,(上三角方程组) （3.2),(n-1) 原方程组化为,以上为消元过程。,(n) 回代求解公式,（3.3) 是回代过程。,（3.3),例1：在MATLAB上，用Gauss列主元消去法求解方程组：,clear a = -0.04 0.04 0.12 3; 0.56 -1.56 0.32 1; -0.24 1.24 -0.28 0 %先

21、定义增广矩阵 x = 0,0,0; %先将解设为零向量，后面重新赋值 tempo = a(2,:); a(2,:) = a(1,:); a(1,:)=tempo; a %第一次选主元（第一行与第二行交换）,a(2,:) = a(2,:) - a(1,:)*a(2,1)/a(1,1); %将第一个对角元下面的数字消为0 a(3,:) = a(3,:) - a(1,:)*a(3,1)/a(1,1); a,a = 0.5600 -1.5600 0.3200 1.0000 0 -0.0714 0.1429 3.0714 0.0000 0.5714 -0.1429 0.4286,tempo = a(3,

22、:); a(3,:) = a(2,:); a(2,:)=tempo;a %第二次选主元 a(3,:) = a(3,:) - a(2,:)*a(3,2)/a(2,2); a %第二次将第二个对角元下的数字消为0,x(3) = a(3,4)/a(3,3); %消去完成，现在开始回代 x(2) = (a(2,4) - a(2,3)*x(3)/a(2,2); x(1) = (a(1,4) - a(1,2:3)*x(2:3)/a(1,1); x,运行得方程组的解为：,a = 0.5600 -1.5600 0.3200 1.0000 0.0000 0.5714 -0.1429 0.4286 0.0000

23、0 0.1250 3.1250,x = 7.0000 7.0000 25.0000,也可直接用 x=Ab,说明： (1)也可采用无回代的列主元消去法(叫Gauss-Jordan消去法)，该法同时消去对角元上下的元素，且仍旧需要选主元，但比有回代的列主元消去法的乘除运算次数多。GaussJordan消去法的优点之一是用它来算逆矩阵的算法非常容易解释。 (2)有回代的列主元消去法所进行的乘除运算次数为 ，计算量很小。,例2：用GaussJordan消去法求解上例中的矩阵 的逆矩阵。,clear A = -0.04 0.04 0.12 ; 0.56 -1.56 0.32 ; -0.24 1.24 -

24、0.28 a=A,eye(3);a tempo = a(2,:); a(2,:) = a(1,:); a(1,:)=tempo; a%第一次选主元，并交换第一行和第二行 a(1,:) = a(1,:)/a(1,1) %将主元标准化 for i=2:3; a(i,:)=a(i,:) - a(i,1)*a(1,:); end; %将主元下的元素消为零 a,tempo = a(3,:); a(3,:) = a(2,:); a(2,:)=tempo; a %第二次选主元，交换第二行和第三行 a(2,:)=a(2,:)/a(2,2); a % 第二个主元标准化,a = 0.5600 -1.5600 0.

25、3200 1.0000 0 -0.0714 0.1429 3.0714 0.0000 0.5714 -0.1429 0.4286,for i=1:3 if i=2, a(i,:)=a(i,:)-a(i,2)*a(2,:); end end % a(2,2)的上下元素消为0 a,a = 1.0000 0 -0.1250 0 3.8750 4.8750 0 1.0000 -0.2500 0 0.7500 1.7500 0 0 0.1250 1.0000 0.1250 0.1250,a(3,:)=a(3,:)/a(3,3) for i=1:3; if i=3, a(i,:)=a(i,:)-a(i,3

26、)*a(3,:); end; end; a A_inv = a(:,4:6) A*A_inv,结果为：A_inv = 1.0000 4.0000 5.0000 2.0000 1.0000 2.0000 8.0000 1.0000 1.0000,也可用rref命令来求,三、 Gauss 全主元消去法： 优点-计算结果更可靠； 缺点-挑主元花时间更多， 次序有变动，程序复杂。,四、Gauss消元法的应用,1、求行列式：det(A),2、求逆矩阵： inv（A）或 A(-1) rref(A,eye(size(A) 将(A,E)化为行最简形，其实就是Gauss-Jordon消去法,（3）求解方程组Ax

27、=b,求逆矩阵的思想: x=inv(A)*b或x=A(-1)*b,左除法（原理上是运用高斯消元法求解,但Matlab在实际执行过程中是通过LU分解法进行的）：x=Ab,符号矩阵法（此法最接近精确值，但运算速度慢）x=sym(A)sym(b),化为行最简形：C=A,b rref(C),例1：在MATLAB上，用Gauss列主元消去法求解方程组：,A= -0.04 0.04 0.12; 0.56 -1.56 0.32; -0.24 1.24 -0.28; b=3;1;0; x11=inv(A)*b %法1：求逆思想 x12=A(-1)*b %法1：求逆思想 x3=Ab %法2：左除法 x4=sym

28、(A)sym(b) %法3：符号矩阵法 C=A,b;rref(C) %法4：化为行最简形,定义3.1,叫,的三角（因子）分解，其中 是,是上三角。,下三角,若L为单位下三角阵（对角元全为1），U为上三角阵，则称 A=LU 为Doolittle分解;,若 L 是下三角，U是单位上三角，则称A=LU为Crout分解。,定理3.1 n 阶阵,有唯一Doolittle分解(Crout）,的所有顺序主子式不为0.,三角分解不唯一,为此引入,定义3.2,五利用矩阵三角分解法求解方程组,为什么要讨论三角分解？若在消元法进行前能实现三角分解,， 则,容易回代求解,在Gauss消去法中，选主元改变了行的次序，尽

29、管对于Gauss消去法来说，这种次序的变换无法事先知道，但是这个变化的影响却可以用一个算子P表示，其中P是一个置换矩阵。用P左乘原始矩阵A得到： PAx=Py 或 对 做Gauss消去法不需要选主元，所以对 做LU分解，同样不需要选主元。,实际上，将选主元Gauss消去法里的行交换同样作用于单位矩阵，所得矩阵即为P。,在MATLAB中，LU分解的命令是lu，有两种格式：,l,u,p=lu(A) 其中A是待分解矩阵；l,u,p分别代表L(主对角线元素为1的下三角矩阵）、U（上三角矩阵）和由单位阵变换出的置换矩阵P 满足：PA=LU，即 A=P-1LU,l,u=lu(A) 其中l= P-1L，所以

30、A=lU=P-1LU 用此格式求出来的 l 并不一定是真正的下三角矩阵， 需要换行后才能是真正的下三角矩阵,实现LU分解后， 若采用l,u=lu(A) ，则Ax=b的解为 x=u (l b) 若采用l，u，p=lu(A)，则Ax=b的解为 x=u( l ( p*b),例 利用LU分解法求解方程组,A=2 4 2 6;4 9 6 15;2 6 9 18;6 15 18 40%先录入系数矩阵 b=9 23 22 47 L,U=lu(A) %对A矩阵做LU分解 y=Lb %求解方程组Ly=bx=Uy %求解方程组Ux=y得到方程组的最终解,故方程组的最终解为: x =(0.5000,2.0000,3

31、.0000,-1.0000),或 A=2 4 2 6;4 9 6 15;2 6 9 18;6 15 18 40%先录入系数矩阵 b=9 23 22 47 L,U,P=lu(A) %对A矩阵做LU分解 y=LP*b %求解方程组Ly=Pbx=Uy %求解方程组Ux=y得到方程组的最终解,X=QR：X为方阵，Q为正交矩阵， R为上三角矩阵,六、利用矩阵QR分解法求解方程组,Q，R=qr(X)：产生一个一个正交矩阵Q和一个上三角矩阵R，满足X=QR,Q，R，E=qr(X)：产生一个一个正交矩阵Q、一个上三角矩阵R、置换矩阵E，满足XE=QR,实现QR分解后， 若采用Q，R=qr(X) ，则Ax=b的

32、解为 x=R(Qb) 若采用Q，R，E=qr(X)，则Ax=b的解为 x=E(R(Qb),例：用QR分解求解线性方程组。 A=3，1，-4，1;1，-3，0，2;0，2，1，-1; 1，6，- 1，-3; b=12，-6，4，0; Q，R=qr(A); x=R(Qb) 或采用QR分解的第2种格式，命令如下： Q，R，E=qr(A); x=E*(R(Qb) 两种得到结果都是 x = -16.4444 20.6667 -1.1111 36.2222,七、利用矩阵Cholesky分解法求解方程组,若矩阵X是对称正定的， X=RR, R为上三角矩阵,R=chol(X)：产生一个上三角阵R，使RR=X

33、若X为非对称正定，则输出一个出错信息,R，p=chol(X)：这个命令格式将不输出出错信息。当X为对称正定的，则p=0，R与上述格式得到的结果相同；否则p为一个正整数。 如果X为满秩矩阵，则R为一个阶数为q=p-1的上三角阵，且满足 RR=X(1:q，1:q),实现Cholesky分解后，线性方程组Ax=b变成RRx=b，所以x=R(Rb)。,【例】用Cholesky分解求解线性方程组。 A=3，1，-4，1;1，-3，0，2;0，2，1，-1;1，6，-1，-3; b=12，-6，4，0; R=chol(A) 结果： ? Error using = chol Matrix must be p

34、ositive definite 命令执行时，出现错误信息，说明A为非正定矩阵。,八、解三对角方程组追赶法,给定方程组,按行严格对角占优,（三对角方程组）,其求解算法是Gauss消去法的一种变形，称为三对角法（追赶法）。解法如下：,1、由第一个方程,令,2、将其代入第二个方程 ，得：,再令,3、将其代入第三个方程,得,再令,以此类推，,令,将,代入最后一个方程,令,以上过程称为 “追”；,总结“追”的算法：,Step1：（初始化变量）,Step2：,下面开始“赶”：,Step3：先求,Step4：再求其他的,三对角方程组的求解程序 tri_diag.m 如下：,% tri_diag(a,b,c

35、,d,n) solves a tridiagonal equation. function x = tri_diag(a,b,c,d,n) for i=2:n r=a(i)/b(i-1); d(i)=d(i)-r*d(i-1); b(i)=b(i)-r*c(i-1); end d(n)=d(n)/b(n); for i=n-1:-1:1 d(i)=(d(i)-c(i)*d(i+1)/b(i); end x=d;,迭代法适用于解大型稀疏方程组(万阶以上的方程组,系数矩阵中零元素占很大比例,而非零元按某种模式分布) 背景: 电路分析、边值问题的数值解和数学物理方程 问题: (1)如何构造迭代格式？

36、 (2)迭代格式是否收敛？ (3)收敛速度如何？ (4)如何进行误差估计？,3.2 解线性方程组的迭代法,一、稀疏矩阵存储方式的产生与转化,1、由完全存储方式 转为 稀疏存储方式命令：,B=sparse(A) 将矩阵A转化为稀疏存储方式的矩阵B sparse(m,n) 生成一个mn的所有元素都是0的稀疏 矩阵 sparse(u,v,S) u,v,S是3个等长的向量。 S是要建立的稀疏矩阵的非0元素； u(i)、v(i)分别是S(i)的行标和列标； 该函数生成一个max(u)行、max(v)列并以S为稀疏元素的稀疏矩阵, u=1,1,4; v=1,2,1; s=1,3,7; sparse(u,v

37、,s),结果：ans = (1,1) 1 (4,1) 7 (1,2) 3,U,V,S=find(A) 返回矩阵A中非0元素的下标和元素，这里产生的U、V、S可作sparse(u,v,s)的参数 full(A) 返回稀疏存储矩阵A对应的完全存储方式矩阵,相关操作的函数：,2、 产生一个稀疏矩阵 把要建立的稀疏矩阵的非0元素及其所在行和列的位置表示出来后由MATLAB自己产生其稀疏存储方式，这需要使用spconvert函数。调用格式为：,B=spconvert(A) 其中A为一个m3或m4的矩阵，每行表示一个非0元素，m是非0元素的个数，A每个元素的意义是：(i,1) 第i个非0元素所在的行。(i

38、,2) 第i个非0元素所在的列。(i,3) 第i个非0元素值的实部。(i,4) 第i个非0元素值的虚部，若矩阵的全部元素都是实数，则无须第四列。, A=1,1,4; 1,2,1; 1,3,7; spconvert(A),结果： ans = (1,1) 4 (1,2) 1 (1,3) 7,3、单位稀疏矩阵的产生 eye 产生一个完全存储方式的单位矩阵 speye(m,n) 产生一个mn的稀疏存储单位矩阵, speye(2,3) ans = (1,1) 1 (2,2) 1,二、简单迭代法 1.迭代法建立. 考虑 Ax=b,(矩阵B不唯一),对应写出,产生向量序列,若收敛,记,则于(1)两端取极限有

39、：,上式说明：,是解向量 ，从而当k充分大时,注意: 迭代阵B不唯一，B的选取影响收敛性。,解向量,(1)叫简单迭代法,B叫迭代矩阵。,2.收敛性. 定义 称,为矩阵B的谱半径。,定理,定理 对于简单迭代法，若迭代矩阵,设有方程组( 其中 ) Ax = b，即,(2),作等价变形,(3),-Jacobi迭代法,(k=0,1,2,),(4),法1：,法2：,Ax = b,A=D+(L+U),Dx = - (L+U ) x + b,= x = - D-1(L+U)x + D-1 b,= x = BJ x+ f,= 迭代公式：x(k+1) = BJ x(k)+ f , (k=0,1,2),BJ= -

40、D-1（L+U）, f =D-1b,编写实现Jacobo迭代法的函数jacobi.m如下:,function s=jacobi(a,b,x0,err) % jacobi迭代法求解线性方程组，a为系数矩阵，b为%ax=b右边的%矩阵b,x0为迭代初值，err为迭代误差 if nargin=3 err=1.0e-6; elseif nargin3 error return end D=diag(diag(a);%构造对角阵D D1=inv(D);%求对角阵D的逆矩阵 L=tril(a,-1);%构造严格下三角阵 U=triu(a,1);%构造严格上三角阵 B=-D1*(L+U);%求出迭代矩阵,f

41、=D1*b; s=B*x0+f; n=1;%n为迭代次数 while norm(s-x0)=err n=n+1 x0=s; s=B*x0+f; s end n,A=9 -1 -1;-1 10 -1;-1 -1 15; b=7;8;13; x0=0;0;0; err=0.00005; s=jacobi(A,b,x0,err),结果：n 1 0.7778 0.8000 0.8667 2 0.9630 0.9644 0.9719 3 0.9929 0.9935 0.9952 4 0.9987 0.9988 0.9991 5 0.9998 0.9998 0.9998 6 1.0000 1.0000 1

42、.0000 7 1.0000 1.0000 1.0000,例,三、 Gauss-Seidel 迭代法,对于,Ax = b,Jacobi迭代公式为,(k=0,1,2,),(4),与（4）对应的迭代公式为：, Gauss-Seidel迭代法,(5),法1：,法2：,Ax = b,A=D+L+U,(D+L) x = - U x + b,= 迭代公式：x(k+1) = BG-S x(k)+ f , (k=0,1,2),编写实现Seidel迭代法的函数seidel.m如下:,function s=seidel(a,b,x0,err) % seidel迭代法求解线性方程组，a为系数矩阵，b为%ax=b右边

43、的%矩阵b,x0为迭代初值，err为迭代误差 if nargin=3 err=1.0e-6; elseif nargin3 error return end D=diag(diag(a);%构造对角阵D L=tril(a,-1);%构造严格下三角阵 U=triu(a,1);%构造严格上三角阵 C=inv(D+L);,B=-C*U; f=C*b; n=1%n为迭代次数 s=B*x0+f while norm(s-x0)=err n=n+1 x0=s; s=B*x0+f; s end n,A=9 -1 -1;-1 10 -1;-1 -1 15; b=7;8;13; x0=0;0;0; err=0.00005; s=seidel(A,b,x0,err),结果：n 1 0.7778 0.8778 0.