第三章 MATLAB线性方程组.ppt

单击此处编辑母版标题样式,*,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,第三章 线性代数方程组及矩阵特征值,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,条对角线以下的元素,上三角矩阵,triu(A,),提取矩阵,A,的上三角矩阵,triu(A,k,),提取矩阵,A,的第,k,条对角线以下的元素,二、,用于专门学科中的矩阵,(1),魔方矩阵,魔方矩阵是每行、每列及两条对角线上的元素和都相等的矩阵。对于,n,阶魔方阵，其元素由,1,2,3,n,2,共,n,2,个整数组成,.,magic(n,),生成一个,n,阶魔方阵,各行各列及两条 对角线和为,(1+2+3+.+n,2,)/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,。命令如下：,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,矩阵的每个元素,h,ij,=1/(i+j-1),hilb(n,),生成,n,阶希尔伯特矩阵,invhilb(n,),专求,n,阶的希尔伯特矩阵的逆矩阵,注意,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 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,，,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=,x,1,x,2,T,1,-1,1,|X|,1,=1,1,1,-1,-1,|X|,2,=1,-1,1,1,-1,-1,|X|,=,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,的,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,矩阵病态特征,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 norm2=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.,Results 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.,Results 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(),会给出正确的结论,例,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.44424930131968 -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)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,-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,),例：,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,是一个,n,n,矩阵，,为,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,的列向量。,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,（,超定方程），多用于曲线拟合,。,解线性方程组的一般函数文件如下：,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,%,方程组不相容,（无解）,，给出最小二乘解,disp,(,方程组的最小二乘法解是：,);,x=Ab;,end,else,%,齐次方程组,if,rank(A,)=n%,列满秩,x=zero(n,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=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%,先定义增广矩阵,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,:);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.1250 3.1250,x=,7.0000,7.0000,25.0000,也可直接用,x=Ab,说明：,(1),也可采用无回代的列主元消去法,(,叫,Gauss-Jordan,消去法,),，该法同时消去对角元上下的元素，且仍旧需要选主元，但比有回代的列主元消去法的乘除运算次数多。,Gauss,Jordan,消去法的优点之一是用它来算逆矩阵的算法非常容易解释。,(2),有回代的列主元消去法所进行的乘除运算次数为 ，计算量很小。,例,2,：用,Gauss,Jordan,消去法求解上例中的矩阵,的逆矩阵。,clear,A=-0.04 0.04 0.12;0.56-1.56 0.32;-0.24 1.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.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)*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=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(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,消去法中，选主元改变了行的次序，尽管对于,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,-1,LU,l,u,=,lu,(,A,),其中,l,=,P,-1,L,，所以,A,=,lU,=P,-1,LU,用此格式求出来的,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.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=,Pb,x=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,的解为,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=R,R,R,为上三角矩阵,R=,chol(X,),：产生一个上三角阵,R,，使,R,R,=,X,若,X,为非对称正定，则输出一个出错信息,R,，,p=,chol(X,),：这个命令格式将不输出出错信息。当,X,为对称正定的，则,p=0,，,R,与上述格式得到的结果相同；否则,p,为一个正整数。,如果,X,为满秩矩阵，则,R,为一个阶数为,q=p-1,的上三角阵，且满足,RR=X(1:q,，,1:q),实现,Cholesky,分解后，线性方程组,Ax=b,变成,RR,x,=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,、将其代入第三个方程,得,再令,以此类推，,令,将,代入最后一个方程,令,以上过程称为“追”；,总结“追”的算法：,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);,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),如何构造迭代格式？,(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,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,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),两端取极限有：,上式说明：,是解向量,，,从而当,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,L,U,=,x=B,J,x+f,=,迭代公式：,x,(k+1),=B,J,x,(k,),+f,(k=0,1,2),B,J,=-D,-,1,（,L+U,）,f,=D,-1,b,D,编写实现,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,nargin,=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.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),=B,G-S,x,(k,),+f,(k=0,1,2),L,U,D,编写实现,Seidel,迭代法的函数,seidel.m,如下,:,function s=seidel(a,b,x0,err),%,seidel,迭代法求解线性方程组，,a,为系数矩阵，,b,为,%ax=b,右边的,%,矩阵,b,x0,为迭代初值，,err,为迭代误差,if,nargin,=3,err=1.0e-6;,elseif,nargin,=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.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,方法对任意 收敛。,法,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,nargin,=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;,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.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