第三章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) 提取矩阵
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 -
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 第三 MATLAB 线性方程组
链接地址:https://www.31doc.com/p-2094887.html