MATLAB线性方程组 (2).ppt

上传人:牧羊曲112 文档编号:6512088 上传时间:2023-11-08 格式:PPT 页数:108 大小:967.50KB
返回 下载 相关 举报
MATLAB线性方程组 (2).ppt_第1页
第1页 / 共108页
MATLAB线性方程组 (2).ppt_第2页
第2页 / 共108页
MATLAB线性方程组 (2).ppt_第3页
第3页 / 共108页
MATLAB线性方程组 (2).ppt_第4页
第4页 / 共108页
MATLAB线性方程组 (2).ppt_第5页
第5页 / 共108页
点击查看更多>>
资源描述

《MATLAB线性方程组 (2).ppt》由会员分享,可在线阅读,更多相关《MATLAB线性方程组 (2).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*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.814749767106

13、56 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.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,1;1,1,3,1;B=pinv(A)例 求矩阵A的伪逆

14、 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=

15、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的根,八、求矩阵A的共轭矩阵 conj(A),eig(A)求矩阵A的特征

16、值 常用的调用格式有: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=

17、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(方程组的最小二乘法

18、解是:);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)分别显示其求解结果。,求解线性方程组(m=n),用克莱姆法则,理论上

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

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

21、,:)=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.0

22、000,也可直接用 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=a(2,:);a(2,:)=a(1,:);a(1,:)=

23、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,:

24、)=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 2.0000 1.0000 2.0000 8.0000 1

25、.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

26、(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 是下三角

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

28、法里的行交换同样作用于单位矩阵,所得矩阵即为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=b的解为 x=u(l(p*b),例 利用LU分解法求解方

29、程组,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=y得到方程组的最终解,X=QR:X为方阵,Q为正交矩阵,R

30、为上三角矩阵,六、利用矩阵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

31、(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=b变成RRx=b,所以x=R(Rb)。,【例】用Chol

32、esky分解求解线性方程组。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、将其代入第三个方程,得,再令,以此类推,,令,将,代入最后一个方程,令,以上过程称为“

33、追”;,总结“追”的算法:,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)=(d(i)-c(i)*d(i+1)/b(i);endx=

34、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)

35、的行标和列标;该函数生成一个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)其

36、中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=

37、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,=迭代公式:x(k+1)=BJ x(k)+f,(k=0,1,2

38、),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);%构造严格上三角阵B=-D1*(L+U);%求出迭代矩阵,f=D1*b;

39、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 1.0000 1.0000 7 1.0000 1.0000

40、 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为迭代初值,err为迭代误差if nargin=3 err=1.0

41、e-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 0.7778 0.8778 0.9770 2 0.9839

42、 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迭代法的函

43、数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*x0+fwhile norm(s-x0)=err n=n

44、+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.3 不可解问题,线性方程组并不总是数值可解的,考虑如下三

45、个方程组。,在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,

展开阅读全文
相关资源
猜你喜欢
相关搜索
资源标签

当前位置:首页 > 生活休闲 > 在线阅读


备案号:宁ICP备20000045号-2

经营许可证:宁B2-20210002

宁公网安备 64010402000987号