《04线性代数中的数值计算.ppt》由会员分享,可在线阅读,更多相关《04线性代数中的数值计算.ppt(46页珍藏版)》请在三一办公上搜索。
1、第4章线性代数中的数值计算,黑龙江大学 电子工程学院,【本章学习目标】,掌握生成特殊矩阵的方法。掌握矩阵分析的方法。掌握求解线性方程组的各种方法。了解矩阵的稀疏存储方式。,目录,4.1 特殊矩阵的生成4.2 矩阵分析4.3 线性方程组求解4.4 矩 阵 分 解4.5 超越函数运算4.6 稀疏矩阵的处理,4.1 特殊矩阵的生成,4.1.1 通用的特殊矩阵zeros函数:产生全0矩阵,即零矩阵。ones函数:产生全1矩阵,即幺矩阵。eye函数:产生单位矩阵,即对角线上的元素为1、其余元素为0的矩阵。rand函数:产生01均匀分布的随机矩阵。randn函数:产生均值为0、方差为1的标准正态分布随机矩
2、阵。这几个函数的调用格式相似,如果这个函数的参数只是一个,那么MATLAB将会创建一个方阵,行数和列数均为这个参数;如果这个函数的参数有两个,那么第一个参数代表行数,第二个参数代表列数。下面以产生零矩阵的zeros函数为例进行说明。zeros函数的调用格式如下。zeros(m):产生mm零矩阵。zeros(m,n):产生mn零矩阵。当m=n时,等同于zeros(m)。zeros(size(A):产生与矩阵A同样大小的零矩阵。,【例4.1】分别建立33、32和与矩阵A同样大小的零矩阵。(1)建立一个33的零矩阵。zeros(3)ans=0 0 0 0 0 0 0 0 0(2)建立一个32的零矩阵
3、。zeros(3,2)(3)设A为23矩阵,则可以用zeros(size(A)建立一个与矩阵A同样大小的零矩阵。A=1 2 3;4 5 6;%产生一个23阶矩阵Azeros(size(A)%产生一个与矩阵A同样大小的零矩阵,【例4.2】建立随机矩阵:(1)在区间10,30内均匀分布的4阶随机矩阵。(2)均值为0.6、方差为0.1的4阶正态分布随机矩阵。产生(0,1)区间均匀分布随机矩阵使用rand函数,假设得到了一组满足(0,1)区间均匀分布的随机数xi,则若想得到在任意a,b区间上均匀分布的随机数,只需用yi=a+(ba)xi计算即可。产生均值为0、方差为1的标准正态分布随机矩阵使用rand
4、n函数,假设已经得到了一组标准正态分布随机数xi,如果想更一般地得到均值为、方差为2的随机数,可用yi=+xi计算出来。针对本例,命令如下:a=10;b=30;x=a+(b-a)*rand(4)y=0.6+sqrt(0.1)*randn(4),4.1.2 面向特定应用的特殊矩阵1魔方矩阵魔方矩阵有一个有趣的性质,其每行、每列及两条对角线上的元素和都相等。对于n阶魔方阵,其元素由1,2,3,n2共n2个整数组成。MATLAB提供了求魔方矩阵的函数magic(n),其功能是生成一个n阶魔方阵。【例4.3】将101125等25个数填入一个5行5列的表格中,使其每行、每列及对角线的和均为565。一个5
5、阶魔方矩阵的每行、每列及对角线的和均为65,对其每个元素都加100后这些和变为565。完成其功能的命令如下:M=100+magic(5)M=117 124 101 108 115 123 105 107 114 116 104 106 113 120 122 110 112 119 121 103 111 118 125 102 109,2范得蒙矩阵范得蒙(Vandermonde)矩阵的最后一列全为1,倒数第二列为一个指定的向量,其他各列是其后列与倒数第二列的点乘积。可以用一个指定向量生成一个范得蒙矩阵。在MATLAB中,函数vander(V)生成以向量V为基础向量的范得蒙矩阵。A=vande
6、r(1:4)A=1 1 1 1 8 4 2 1 27 9 3 1 64 16 4 1,3希尔伯特矩阵希尔伯特(Hilbert)矩阵是一种数学变换矩阵,它的每个元素hij=1/(i+j1)。在MATLAB中,生成希尔伯特矩阵的函数是hilb(n)。专门求希尔伯特矩阵的逆的函数invhilb(n)format rat%以有理形式输出H=hilb(4)H=1 1/2 1/3 1/4 1/2 1/3 1/4 1/5 1/3 1/4 1/5 1/6 1/4 1/5 1/6 1/7 H=invhilb(4)format short%恢复默认输出格式,4托普利兹矩阵托普利兹(Toeplitz)矩阵除第一行第
7、一列外,其他每个元素都与左上角的元素相同。生成托普利兹矩阵的函数是toeplitz(x,y),它生成一个以x为第一列、y为第一行的托普利兹矩阵。这里x、y均为向量,两者不必等长。toeplitz(x)用向量x生成一个对称的托普利兹矩阵。T2=toeplitz(1:4)T2=1 2 3 4 2 1 2 3 3 2 1 2 4 3 2 1,5伴随矩阵生成伴随矩阵的函数是compan(p),其中p是一个多项式的系数向量,高次幂系数排在前,低次幂排在后。例如,为了求多项式的x37x+6的伴随矩阵,可使用如下命令:p=1,0,-7,6;compan(p)ans=0 7-6 1 0 0 0 1 0,6帕斯
8、卡矩阵我们知道,二次项(x+y)n展开后的系数随n的增大组成一个三角形表,称为杨辉三角形。由杨辉三角形表组成的矩阵称为帕斯卡(Pascal)矩阵,它的元素p1j=1,p i1=1,p ij=p i1,j1+p i1,j(i1,j1)。函数pascal(n)生成一个n阶帕斯卡矩阵。【例4.5】求(x+y)4的展开式。在MATLAB命令窗口,输入命令:pascal(5)ans=1 1 1 1 1 1 2 3 4 5 1 3 6 10 15 1 4 10 20 35 1 5 15 35 70矩阵次对角线上的元素1,4,6,4,1即为展开式的系数,即(x+y)4=x4+4x3y+6x2y2+4xy3+
9、y4,4.2 矩阵分析,4.2.1 矩阵结构变换1对角阵只有对角线上有非0元素的矩阵称为对角矩阵,对角线上的元素都为1的对角矩阵称为单位矩阵。(1)提取矩阵的对角线元素设A为mn矩阵,函数diag(A)用于提取矩阵A主对角线元素,产生一个具有min(m,n)个元素的列向量。例如:A=1,2,3;4,5,6A=1 2 3 4 5 6D=diag(A)D=1 5diag(A,k)提取第k条对角线的元素。主对角线为第0条对角线;与主对角线平行,往上为第1条,第2条,第n条对角线,往下为第1条,第2条,第n条对角线。,(2)构造对角矩阵设V为具有m个元素的向量,diag(V,k)的功能是产生一个nn(
10、n=m+|k|)对角阵,其第k条对角线的元素即为向量V的元素。例如:diag(1:3,-1)ans=0 0 0 0 1 0 0 0 0 2 0 0 0 0 3 0省略k时,相当于k为0,其主对角线元素即为向量V的元素。,【例4.6】先建立55矩阵A,然后将A的第一行元素乘以1,第二行乘以2,第五行乘以5。用一个对角矩阵左乘一个矩阵时,相当于用对角阵的第一个元素乘以该矩阵的第一行,用对角阵的第二个元素乘以该矩阵的第二行依此类推,因此,只需按要求构造一个对角矩阵D,并用D左乘A即可。命令如下:A=1:5;2:6;3:7;4:8;5:9D=diag(1:5);D*A%用D左乘A,对A的每行乘以一个指
11、定常数,2三角阵三角阵又进一步分为上三角阵和下三角阵。所谓上三角阵,即矩阵的对角线以下的元素全为0的一种矩阵,而下三角阵则是对角线以上的元素全为0的一种矩阵。与矩阵A对应的上三角阵B是与A具有相同的行数和列数的一个矩阵,并且B的对角线以上(含对角线)的元素和A对应相等,而对角线以下的元素等于0。求矩阵A的上三角阵的MATLAB函数是triu(A)。triu(A,k),其功能是求矩阵A的第k条对角线以上的元素。在MATLAB中,提取矩阵A的下三角矩阵的函数是tril(A)和tril(A,k),3矩阵的转置所谓转置,即把源矩阵的第一行变成目标矩阵第一列,第二行变成第二列依此类推。显然,一个m行n列
12、的矩阵经过转置运算后,变成一个n行m列的矩阵。MATLAB中,转置运算符是单撇号()。A=1:3;1:3;1:3B=A,4矩阵的旋转在MATLAB中,可以很方便地以90为单位对矩阵A按逆时针方向进行旋转。利用函数rot90(A,k)将矩阵A旋转90的k倍,当k为负整数时,对矩阵A按顺时针方向进行旋转;当k为1时可省略。例如,将A按逆时针方向旋转90,命令如下:A=9,37,38;-2,31,8;0,84,5;B=rot90(A)B=38 8 5 37 31 84 9-2 0rot90(A,4)ans=9 37 38-2 31 8 0 84 5,5矩阵的翻转矩阵的翻转分左右翻转和上下翻转。对矩阵
13、实施左右翻转是将原矩阵的第一列和最后一列调换,第二列和倒数第二列调换依此类推。对矩阵A实施左右翻转的函数是fliplr(A)。对矩阵A实施上下翻转的函数是flipud(A)。,4.2.2 矩阵求值1方阵的行列式值把一个方阵看作一个行列式,并对其按行列式的规则求值,这个值就称为矩阵所对应的行列式的值。在MATLAB中,求方阵A所对应的行列式的值的函数是det(A)。例如:A=1:3;2:-1:0;12,5,9A=1 2 3 2 1 0 12 5 9B=det(A)B=-33,2矩阵的秩与迹(1)矩阵的秩矩阵线性无关的行数和列数称为矩阵的秩。rank(A)A=1,2,3;2,5,6;3,2,5;r
14、=rank(A)r=3(2)矩阵的迹矩阵的迹即矩阵的对角线元素之和。trace(A)。,3向量和矩阵的范数,函数norm用于计算矩阵或向量的范数,norm函数的格式如下。norm(X,1):求向量或矩阵X的1范数。norm(X)、norm(X,2):求向量或矩阵X的2范数。norm(X,inf):求向量或矩阵X的范数。,4矩阵的条件数矩阵A的条件数等于A的范数与A的逆矩阵的范数的乘积,即。这样定义的条件数总是大于1的。条件数越接近于1,矩阵的性能越好,反之,矩阵的性能越差。A有3种范数,相应地可定义3种条件数。在MATLAB中,计算A的3种条件数的函数如下。cond(A,1):计算A的1范数下
15、的条件数。cond(A)或cond(A,2):计算A的2范数数下的条件数。cond(A,inf):计算A的范数下的条件数。,4.2.3 矩阵的特征值与特征向量对于n阶方阵A,求数和向量,使得等式A=成立,满足等式的数称为A的特征值,而向量称为A的特征向量。在MATLAB中,计算矩阵A的特征值和特征向量的函数是eig(A),常用的调用格式有如下3种。E=eig(A):求矩阵A的全部特征值,构成向量E。V,D=eig(A):求矩阵A的全部特征值,构成对角阵D,并求A的特征向量构成V的列向量。V,D=eig(A,nobalance):与第2种格式类似,但第2种格式中先对A作相似变换后求矩阵A的特征值
16、和特征向量,而格式3直接求矩阵A的特征值和特征向量。一个矩阵的特征向量有无穷多个,eig函数只找出其中的n个,A的其他特征向量,均可由这n个特征向量的线性组合表示。,例如:A=1,1,0.5;1,1,0.25;0.5,0.25,2;V,D=eig(A)V=0.7212 0.4443 0.5315-0.6863 0.5621 0.4615-0.0937-0.6976 0.7103D=-0.0166 0 0 0 1.4801 0 0 0 2.5365,4.3 线性方程组求解,4.3.1 矩阵求逆及线性代数方程组求解1矩阵求逆若方阵A、B满足等式:A B=B A=I(I为单位阵),称A为B的逆矩阵函
17、数inv(A)用于计算方阵的逆矩阵。A=1-1 1;5-4 3;2 1 1;B=inv(A)B=-1.4000 0.4000 0.2000 0.2000-0.2000 0.4000 2.6000-0.6000 0.2000A*Bans=1.0000 0 0-0.0000 1.0000 0-0.0000 0 1.0000,2利用矩阵求逆方法解线性方程组线性方程组的矩阵表示形式为Ax=b在方程组两边各左乘A-1,有A-1 Ax=A-1 b得到x=A-1 b【例4.9】利用矩阵求逆方法解线性方程组A=1,-2,3;3,-1,5;2,1,5;b=1;2;3;x=inv(A)*bx=-0.3333 0.
18、3333 0.6667,4.3.2 利用左除运算符求解线性方程组对于线性方程组Ax=b,可以利用左除运算符“”求解:x=Ab【例4.10】用左除运算符求解下列相同系数矩阵的两个线性代数方程组的解。,解法1:分别解线性方程组。A=1,-1,1;5,-4,3;2,1,1;b1=2;-3;1;b2=3;4;-5;x=Ab1y=Ab2,解法2:将两个线性方程组连在一起求解。A=1,-1,1;5,-4,3;2,1,1;b=2,3;-3,4;1,-5;xy=Ab,4.4 矩 阵 分 解,矩阵分解是指根据一定的原理用某种算法将一个矩阵分解成若干个矩阵的乘积。常见的矩阵分解有LU分解、QR分解、Cholesk
19、y分解以及奇异分解等。4.4.1 矩阵的LU分解矩阵的LU分解又称Gauss消去分解或三角分解,就是将一个方阵表示为一个行交换下三角矩阵和一个上三角矩阵的乘积形式。MATLAB提供的lu函数用于对矩阵进行LU分解,其调用格式如下。L,U=lu(X):产生一个上三角阵U和一个变换形式的下三角阵L(行交换),使之满足X=LU。注意,这里的矩阵X必须是方阵。L,U,P=lu(X):产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P,使之满足PX=LU。当然矩阵X同样必须是方阵。当使用第1种格式时,矩阵L往往不是一个下三角矩阵,但可以通过行交换成为一个下三角阵。,利用第2种格式对矩阵A进行LU分解:
20、L,U,P=lu(A),clearA=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4;b=13,-9,6,0;ticx2=Ab;%用左除运算求解tocticx1=inv(A)*b;%用求逆运算求解tocticL,U=lu(A);%LU分解x3=U(Lb);%用LU分解求解tocx1=x1x2=x2x3=x3,其中tic和toc两个函数配合使用用于计算程序的执行时间,tic记录当前时间,toc记录或显示使用tic函数以来所花费的时间。运行结果说明,x1、x2、x3的值相同,通过LU分解求值所化运行时间最少。,4.4.2 矩阵的QR分解对矩阵X进行QR分解,就是把X分解为
21、一个正交矩阵Q和一个上三角矩阵R的乘积形式。qr函数可用于对矩阵进行QR分解,其调用格式如下:Q,R=qr(X):产生一个一个正交矩阵Q和一个上三角矩阵R,使之满足X=QR。Q,R,E=qr(X):产生一个一个正交矩阵Q、一个上三角矩阵R以及一个置换矩阵E,使之满足XE=QR。,设对矩阵A进行QR分解。A=2,1,1,4;1,2,-1,2;1,-1,3,3;Q,R=qr(A)Q=-0.8165 0-0.5774-0.4082-0.7071 0.5774-0.4082 0.7071 0.5774R=-2.4495-1.2247-1.6330-5.3072 0-2.1213 2.8284 0.70
22、71 0 0 0.5774 0.5774为验证结果是否正确,输入命令QR=Q*R,4.4.3 矩阵的Cholesky分解如果矩阵X是对称正定的,则Cholesky分解将矩阵X分解成一个下三角矩阵和上三角矩阵的乘积。设上三角矩阵为R,则下三角矩阵为其转置,即X=RR。MATLAB函数chol(X)用于对矩阵X进行Cholesky分解,其调用格式如下。R=chol(X):产生一个上三角阵R,使RR=X。若X为非对称正定,则输出一个出错信息。R,p=chol(X):这个命令格式将不输出出错信息。当X为对称正定的,则p=0,R与上述格式得到的结果相同;否则p为一个正整数。如果X为满秩矩阵,则R为一个阶
23、数为q=p1的上三角阵,且满足RR=X(1:q,1:q)。,设对矩阵A进行Cholesky分解。A=2,1,1;1,2,-1;1,-1,3;R=chol(A)R=1.4142 0.7071 0.7071 0 1.2247-1.2247 0 0 1.0000为验证结果是否正确,输入命令R*Rans=2.0000 1.0000 1.0000 1.0000 2.0000-1.0000 1.0000-1.0000 3.0000,4.5 超越函数运算,1超越函数MATLAB还提供了一些直接作用于矩阵的超越函数,如矩阵平方根函数sqrtm、矩阵指数函数expm、矩阵对数函数logm等,这些函数名都在上述内
24、部函数名之后缀以m,并规定输入参数A必须是方阵。例如:A=4,2;3,6;B=sqrtm(A)B=1.9171 0.4652 0.6978 2.38232通用矩阵函数funmfunm(A,fun)对方阵A计算由fun定义的函数的矩阵函数值。例如,当fun取exp时,funm(A,exp)可以计算矩阵A的指数,与expm(A)的计算结果一样。,例子A=2,-1;1,0;funm(A,exp)ans=5.4366-2.7183 2.7183 0.0000expm(A)ans=5.4366-2.7183 2.7183-0.0000注意:funm函数可以用于exp函数和log函数,但求矩阵的平方根只能
25、用sqrtm。,4.6 稀疏矩阵的处理,稀疏矩阵中具有大量的零元素,而仅含极少量的非零元素。4.6.1 矩阵存储方式1完全存储方式:矩阵的全部元素按列存储。以前讲到的矩阵的存储方式都是按这个方式存储的,此存储方式对稀疏矩阵也适用。例如,不论是mn阶普通的还是稀疏的实矩阵均需要mn个存储单元,而复矩阵还要翻倍。在这种方式下,矩阵中的全部零元素也必须输入。2稀疏存储方式:仅存储矩阵所有的非零元素的值及其位置,即行号和列号。,4.6.2 矩阵的稀疏存储方式1将完全存储方式转化为稀疏存储方式函数A=sparse(S)将矩阵S转化为稀疏存储方式的矩阵A。当矩阵S是稀疏存储方式时,则函数调用相当于A=S。
26、,X=2,0,0,0,0;0,0,0,0,0;0,0,0,5,0;0,1,0,0,-1;0,0,0,0,-5;A=sparse(X)A=(1,1)2(4,2)1(3,4)5(4,5)-1(5,5)-5,sparse函数还有其他一些调用格式。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,v,S=find(A):返回矩阵A中非0元素的下标和元素。这
27、里产生的u、v、S可作为sparse(u,v,S)的参数。full(A):返回和稀疏存储矩阵A对应的完全存储方式矩阵。,2产生稀疏存储矩阵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元素值的虚部,若矩阵的全部元素都是实数,则无须第4列。该函数将A所描述的一个稀疏矩阵转化为一个稀疏存储矩阵。,设命令如下A=2,2,1;3,1,-1;4,3,3;5,3,8;6,6,12;B=spconvert(A)
28、B=(3,1)-1(2,2)1(4,3)3(5,3)8(6,6)12,3带状稀疏存储矩阵,将带状对角线之值构成下列矩阵B,将带状的位置构成向量d:,然后利用spdiags函数产生一个稀疏存储矩阵。B=0,0,0,41,51;11,21,31,42,52;12,22,32,0,0;d=-3,0,3;A=spdiags(B,d,5,6)%产生一个稀疏存储矩阵A,A=spdiags(B,d,m,n)其中,m、n为原带状矩阵的行数与列数。B为rp阶矩阵,这里r=min(m,n),p为原带状矩阵所有非零对角线的条数,矩阵B的第i列即为原带状矩阵的第i条非零对角线。取值方法是:若非零对角线上元素个数等于r
29、,则取全部元素;若非零对角线上元素个数小于r,则应该用零补足到r个元素。补零的原则:当mn,应从该对角线的第1行开始补零或向后补零至末行;当mn,则应从该对角线的第1列开始补零或向后补零至末列。d为具有p个元素的列向量,它的第i个元素为该带状矩阵的第i条对角线的位置k。k的取法:若是主对角线,取k=0,若位于主对角线的下方第s条对角线,取k=s,若位于主对角线的上方第s条对角线,则取k=s。,4单位矩阵的稀疏存储单位矩阵只有对角线元素为1,其他元素都为0,是一种具有稀疏特征的矩阵。函数eye产生一个完全存储方式的单位矩阵。MATLAB还有一个产生稀疏存储方式的单位矩阵的函数,这就是speye。
30、函数speye(m,n)返回一个mn的稀疏存储单位矩阵。例如:s=speye(3,5)s=(1,1)1(2,2)1(3,3)1,4.6.3 稀疏矩阵应用举例稀疏存储矩阵只是矩阵的存储方式不同,它的运算规则与普通矩阵是一样的。当参与运算的对象不全是稀疏存储矩阵时,所得结果一般是完全存储形式。例如:A=0,0,3;0,5,0;0,0,9;B=sparse(A);B*B%两个稀疏矩阵相乘,结果仍为稀疏矩阵ans=(2,2)25(1,3)27(3,3)81rand(3)*B%普通矩阵与随机矩阵相乘,结果为完全存储矩阵ans=0 1.7643 2.4808 0 4.0658 4.5058 0 0.0493 1.9622,