MATLAB线性方程组.ppt

《MATLAB线性方程组.ppt》由会员分享，可在线阅读，更多相关《MATLAB线性方程组.ppt（108页珍藏版）》请在三一办公上搜索。

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)提取矩阵A的第k条对角线以下的元素上三

2、角矩阵 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行5列的表格中，使其每行每列及对角线的和均为565。命令

3、如下：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阶的希尔伯特矩阵的逆矩阵 注意1：高阶Hilbert矩阵一般为病态矩阵，所以直接求逆可能出现错误结论，故应该

4、采用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 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

5、 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，10，10，5，1即为展开式的系数。,三、向量和矩阵的范数 norm(V)或norm(V,2)求向量V（或矩阵V)的2范

6、数 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;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的

7、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)可能出错),例：Hilbert矩阵（非常有名的病态矩阵）：,验证从55到1414的Hilbert矩阵病态特征,clearformat long;for n=5:14 H=hilb

8、(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 norm2=1.455203e-011 n=7 norm1=9.810538e-008 norm2=5.208793e-010 n=8 norm1=3.123671e-006 norm2=4.82218

9、7e-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-003Warning:Matrix is close to singular or badly scaled.Results may be inaccurate.RCOND=2.632766e-017.n=12 norm1=1.176913e+001 norm2=5.576679e-002Warning:Matrix is

10、close to singular or badly scaled.Results may be inaccurate.RCOND=2.339949e-018.n=13 norm1=5.323696e+001 norm2=1.137063e+001Warning:Matrix is close to singular or badly scaled.Results may be inaccurate.RCOND=1.708191e-019.n=14 norm1=1.224232e+004 norm2=2.836298e+002,说明1:对于Hilbert求逆时，不建议用inv,可用 invhi

11、lb直接产生逆矩阵说明2：符号工具箱中也对符号矩阵定义了inv()函数,即使对更高阶的奇异矩阵也可以精确的求解出逆矩阵,例：H=sym(hilb(30);norm(double(H*inv(H)-eye(size(H)结果为ans=0,说明3：对于奇异矩阵用数值方法的inv()函数，会给出错误的结果，但符号工具箱中的inv()会给出正确的结论,例 A=16,2,3,13;5,11,10,8;9,7,6,12;4,14,15,1;det(A)B=inv(A),结果：ans=0Warning:Matrix is close to singular or badly scaled.Results m

12、ay be inaccurate.RCOND=1.306145e-017.B=1.0e+014*,但用符号工具箱的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.invError,(in inverse)singular matrix,2、矩阵的伪逆 pinv(A)若A不是一个方阵，或A是一个非满秩的方阵时，矩阵A没有逆矩阵，但可以找到一个与A的转置矩阵A同型的矩阵B，使得：ABA=A,BAB=B此时称矩阵B为矩阵A的伪逆。例 求A的伪逆，并将结果送BA=3,1,1,1;1,3,1

13、,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,-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;%用

14、方程组的右端向量置换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)例：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的根,八、求矩

15、阵A的共轭矩阵 conj(A),eig(A)求矩阵A的特征值 常用的调用格式有：E=eig(A)求矩阵A的全部特征值，构成向量E。V,D=eig(A)求矩阵A的全部特征值，构成对角阵D，并求A的特征向量构成V的列向量。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

16、(A,nobalance),现求解线性方程组,情形1：m=n（恰定方程组）在MATLAB中的求解命令有：,情形2：mn（超定方程），多用于曲线拟合。,解线性方程组的一般函数文件如下：function 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%方程组不

17、相容（无解），给出最小二乘解 disp(方程组的最小二乘法解是：);x=Ab;endelse%齐次方程组 if rank(A)=n%列满秩 x=zero(n,1)%0解 else%非0解，无穷多个解 disp(方程组有无穷个解，基础解系为x);x=null(A,r);end endreturn,如在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)分别显示其

18、求解结果。,求解线性方程组（m=n),用克莱姆法则，理论上可行，但实际运算时工作量大，耗时。下面研究,、一、Gauss简单消元法（m=n),设 有,用线性代数中的克莱姆法则求解时，工作量非常大，工作量小的方法是 Gauss消元法。,3.1 解线性方程组的直接法,消 元：,以此类推，最后方程组化为：,回 代：,失效，故应选主元,二、列主元素消去法-计算结果可靠,到此原方程组化为,到此原方程组化为,(上三角方程组)（3.2),(n-1)原方程组化为,以上为消元过程。,(n)回代求解公式,（3.3)是回代过程。,（3.3),例1：在MATLAB上，用Gauss列主元消去法求解方程组：,cleara=

19、-0.04 0.04 0.12 3;0.56-1.56 0.32 1;-0.24 1.24-0.28 0%先定义增广矩阵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);%将第一个对角元下面的数字消为0a(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

20、.1429 0.4286,tempo=a(3,:);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 0 0.125

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

22、=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.3200 1.0000 0-0.0714 0.1429 3.0714 0.0000 0.5714-0.1429

23、0.4286,for i=1:3 if i=2,a(i,:)=a(i,:)-a(i,2)*a(2,:);endend%a(2,2)的上下元素消为0a,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)*a(3,:);end;end;aA_inv=a(:,4:6)A*A_inv,结果为：A_inv=1.0000 4.0000 5.0000

24、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=b,求逆矩阵的思想:x=inv(A)*b或x=A(-1)*b,左除法（原理上是运用高斯消元法求解,但Matlab在实际执行过程中是通过LU分解法进行的）：x=Ab,符号矩阵

25、法（此法最接近精确值，但运算速度慢）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(A)sym(b)%法3：符号矩阵法C=A,b;rref(C)%法4：化为行最简形,定义3.1,叫,的三角（因子）分解，其中 是,是上三角。,下三角,若L为单位下三角阵（对角元全为1），U为上三角阵，则

26、称 A=LU 为Doolittle分解;,若 L 是下三角，U是单位上三角，则称A=LU为Crout分解。,定理3.1 n 阶阵,有唯一Doolittle分解(Crout）,的所有顺序主子式不为0.,三角分解不唯一,为此引入,定义3.2,五利用矩阵三角分解法求解方程组,为什么要讨论三角分解？若在消元法进行前能实现三角分解,，则,容易回代求解,在Gauss消去法中，选主元改变了行的次序，尽管对于Gauss消去法来说，这种次序的变换无法事先知道，但是这个变化的影响却可以用一个算子P表示，其中P是一个置换矩阵。用P左乘原始矩阵A得到：PAx=Py 或 对 做Gauss消去法不需要选主元，所以对 做L

27、U分解，同样不需要选主元。,实际上，将选主元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，所以A=lU=P-1LU 用此格式求出来的 l 并不一定是真正的下三角矩阵，需要换行后才能是真正的下三角矩阵,实现LU分解后，若采用l,u=lu(A)，则Ax=b的解为 x=u(l b)若采用l，u，p=lu(A)，则Ax

28、=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 47L,U=lu(A)%对A矩阵做LU分解y=Lb%求解方程组Ly=bx=Uy%求解方程组Ux=y得到方程组的最终解,故方程组的最终解为:x=(0.5000,2.0000,3.0000,-1.0000),或A=2 4 2 6;4 9 6 15;2 6 9 18;6 15 18 40%先录入系数矩阵b=9 23 22 47L,U,P=lu(A)%对A矩阵做LU分解y=LP*b%求解方程组Ly=Pbx=Uy%求解方程组Ux=

29、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的解为 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分解的

30、第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 若X为非对称正定，则输出一个出错信息,R，p=chol(X)：这个命令格式将不输出出错信息。当X为对称正定的，则p=0，R与上述格式得到的结果相同；否则p为一个正整数。如果X为满秩矩阵，则R为一个阶数为q=p-1的上三角阵，且满足 RR=X(1:q，1:q),实现Cholesky分解后，线性方程组Ax

31、=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 positive definite命令执行时，出现错误信息，说明A为非正定矩阵。,八、解三对角方程组追赶法,给定方程组,按行严格对角占优,（三对角方程组）,其求解算法是Gauss消去法的一种变形，称为三对角法（追赶法）。解法如下：,1、由第一个方程,令,2、将其代入第二个方程，得：,再令,3、将其代入第三个方程,得,再令

32、,以此类推，,令,将,代入最后一个方程,令,以上过程称为“追”；,总结“追”的算法：,Step1：（初始化变量）,Step2：,下面开始“赶”：,Step3：先求,Step4：再求其他的,三对角方程组的求解程序 tri_diag.m 如下：,%tri_diag(a,b,c,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);endd(n)=d(n)/b(n);for i=n-1:-1:1 d(i)

33、=(d(i)-c(i)*d(i+1)/b(i);endx=d;,迭代法适用于解大型稀疏方程组(万阶以上的方程组,系数矩阵中零元素占很大比例,而非零元按某种模式分布)背景:电路分析、边值问题的数值解和数学物理方程问题:(1)如何构造迭代格式？(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是要

34、建立的稀疏矩阵的非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,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自己产生其稀疏存储方式，这需要使用spcon

35、vert函数。调用格式为：,B=spconvert(A)其中A为一个m3或m4的矩阵，每行表示一个非0元素，m是非0元素的个数，A每个元素的意义是：(i,1)第i个非0元素所在的行。(i,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

36、(2,2)1,二、简单迭代法 1.迭代法建立.考虑 Ax=b,(矩阵B不唯一),对应写出,产生向量序列,若收敛,记,则于(1)两端取极限有：,上式说明：,是解向量，从而当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,=迭

37、代公式：x(k+1)=BJ x(k)+f,(k=0,1,2),BJ=-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 returnendD=diag(diag(a);%构造对角阵DD1=inv(D);%求对角阵D的逆矩阵L=tril(a,-1);%构造严格下三角阵U=triu(a,1);%构造严格上三

38、角阵B=-D1*(L+U);%求出迭代矩阵,f=D1*b;s=B*x0+f;n=1;%n为迭代次数while norm(s-x0)=err n=n+1 x0=s;s=B*x0+f;sendn,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

39、1.0000 1.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右边的%矩阵b,x0为迭代初

40、值，err为迭代误差if nargin=3 err=1.0e-6;elseif nargin3 error returnendD=diag(diag(a);%构造对角阵DL=tril(a,-1);%构造严格下三角阵U=triu(a,1);%构造严格上三角阵C=inv(D+L);,B=-C*U;f=C*b;n=1%n为迭代次数s=B*x0+fwhile norm(s-x0)=err n=n+1 x0=s;s=B*x0+f;sendn,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

41、0.7778 0.8778 0.9770 2 0.9839 0.9961 0.9987 2 0.9994 0.9998 0.9999 3 1.0000 1.0000 1.0000 4 1.0000 1.0000 1.0000,例,发现：seidel收敛速度比jacobi快,松弛因子，=1 即Seidel方法(5),(6)是一种加权平均。,四、超松弛迭代法(SOR法),法1,SOR方法的收敛性如下（不加证明）：,(1)SOR方法对任意 都收敛的必要条件是：(2)若系数矩阵A对称正定，则 时SOR方法 求解 对任意 收敛；(3)若系数矩阵A按行（或按列）严格对角占优，则 时SOR方法对任意 收敛。

42、,法2：SOR法的矩阵迭代格式：,编写实现SOR迭代法的函数SOR.m如下:,function s=SOR(a,b,x0,w,err)%seidel迭代法求解线性方程组，a为系数矩阵，b为%ax=b右边的%矩阵b,x0为迭代初值，w为松弛因子，%err为迭代误差if nargin=4 err=1.0e-6;elseif nargin4 error returnendD=diag(diag(a);%构造对角阵DL=tril(a,-1);%构造严格下三角阵U=triu(a,1);%构造严格上三角阵C=inv(D+w*L);,B=C*(1-w)*D-w*U;f=w*C*b;n=1%n为迭代次数s=B

43、*x0+fwhile norm(s-x0)=err n=n+1 x0=s;s=B*x0+f;sendn,A=9-1-1;-1 10-1;-1-1 15;b=7;8;13;x0=0;0;0;err=0.00005;w=1.1;s=seidel(A,b,x0,w,err),结果：n 1 0.8556 0.9741 1.0875 2 1.0220 1.0146 0.9939 3 0.9988 0.9977 1.0004 4 0.9999 1.0003 1.0000 5 1.0000 1.0000 1.0000 6 1.0000 1.0000 1.0000,例,发现：SOR收敛速度比jacobi快,3

44、.3 不可解问题,线性方程组并不总是数值可解的，考虑如下三个方程组。,在MATLAB中分别求解如下：,有无穷多个解,其中一个解,方程组不相容,最小二乘解,3.4 病态问题,有许多线性方程组理论上是可解的，但实际计算中由于受到舍入误差的干扰而无法得到精确解，此类问题称为病态问题。,病态矩阵的一个重要标志是条件数,矩阵A的条件数记为Cond(A)，定义为：,条件数总满足：,注：当矩阵是病态时，其条件数一定很大，但并不能直接说明解的误差。,MATLAB中计算条件数的命令是：cond(A),对于病态矩阵，逆矩阵和行列式的计算都会变得不精确。所以具备下列特征的问题可认为是病态的：,下面以一个4阶hilbert矩阵为例，看看病态矩阵对线性方程组解的影响：,ans=-2.4000 27.0000-64.8000 42.0000 ans=524.0568 1.5514e+004 4.7661e+005 1.4951e+007 4.7537e+008 1.5258e+010,A=hilb(4);b=1 2 1.41 2;b1=1 2 1.42 2;Ab-Ab1for k=3:8 H=hilb(k);cond(H)end,