第九讲_线性代数中的数值计算问题

第九讲matlab求解线性代数【引例】求下列三阶线性代数方程组的近似解5426255452321321321xxxxxxxxxMATLAB程序为：A=[2-54;15-2;-124]；b=[5;6;5]；x=A\b在MATLAB命令窗口，先输入下列命令构造系数矩阵A和右端向量b：A=[2-54;15-2;-124]A=2-5415-2-124b=[5;6;5]b=565然后只需输入命令x=A\b即可求得解x：x=A\bx=2.76741.18601.3488一、特殊矩阵的实现1.零矩阵:所有元素值为零的矩阵称为零矩阵。零矩阵可以用zeros函数实现。zeros是MATLAB内部函数，使用格式如下：zeros(m)：产生m阶零矩阵；zeros(m,n)：产生mxn阶零矩阵，当m=n时等同于zeros(m)；zeros(size(A))：产生与矩阵A同样大小的零矩阵。一、特殊矩阵的实现常见的特殊矩阵有零矩阵、幺矩阵、单位矩阵、三角形矩阵等，这类特殊矩阵在线性代数中具有通用性；还有一类特殊矩阵在专门学科中有用，如有名的希尔伯特(Hilbert)矩阵、范德蒙(Vandermonde)矩阵等。2.幺矩阵:所有元素值为1的矩阵称为幺矩阵。幺矩阵可以用ones函数实现。它的调用格式与zeros函数一样。【例1】试用ones分别建立32阶幺矩阵、和与前例矩阵A同样大小的幺矩阵。用ones(3,2)建立一个32阶幺阵：ones(3,2)%一个32阶幺阵ans=111111一、特殊矩阵的实现3.单位矩阵:主对角线的元素值为1、其余元素值为0的矩阵称为单位矩阵。它可以用MATLAB内部函数eye建立，使用格式与zeros相同。4.数量矩阵:主对角线的元素值为一常数d、其余元素值为0的矩阵称为数量矩阵。显然，当d=1时，即为单位矩阵，故数量矩阵可以用eye(m)*d或eye(m,n)*d建立。5.对角阵:对角线的元素值为常数、其余元素值为0的矩阵称为对角阵。我们可以通过MATLAB内部函数diag，利用一个向量构成对角阵；或从矩阵中提取某对角线构成一个向量。使用格式为diag(V)和diag(V,k)两种。6.用一个向量V构成一个对角阵设V为具有m个元素的向量，diag(V)将产生一个m阶对角阵，其主对角线的元素值即为向量的元素值；diag(V,k)将产生一个nxn(n=m+|k|，k为一整数)阶对角阵，其第k条对角线的元素值即为向量的元素值。注意：当k＞0，则该对角线位于主对角线的上方第k条；当k＜0，该对角线位于主对角线的下方第|k|条；当k=0，则等同于diag(V)。用diag建立的对角阵必是方阵。一、特殊矩阵的实现【例2】已知向量v，试建立以向量v作为主对角线的对角阵A；建立分别以向量v作为主对角线两侧的对角线的对角阵B和C。MATLAB程序如下：v=[1;2;3];%建立一个已知的向量AA=diag(v)A=100020003B=diag(v,1)B=0100002000030000C=diag(v,-1)C=0000100002000030%按各种对角线情况构成相应的对角阵A、B和C7.从矩阵中提取某对角线我们也可以用diag从矩阵中提取某对角线构成一个向量。设A为mn阶矩阵，diag(A)将从矩阵A中提取其主对角线产生一个具有min(m,n)个元素的向量。diag(A,k)的功能是：当k＞0，则将从矩阵A中提取位于主对角线的上方第k条对角线构成一个具有n-k个元素的向量；当k＜0，则将从矩阵A中提取位于主对角线的下方第|k|条对角线构成一个具有m+k个元素的向量；当k=0，则等同于diag(A)。【例3】已知矩阵A，试从矩阵A分别提取主对角线及它两侧的对角线构成向量B、C和D。MATLAB程序如下：A=[123;456];%建立一个已知的23阶矩阵A%按各种对角线情况构成向量B、C和DB=diag(A)B=15C=diag(A,1)C=26D=diag(A,-1)D=48.上三角阵：使用格式为triu(A)、triu(A,k)设A为mn阶矩阵，triu(A)将从矩阵A中提取主对角线之上的上三角部分构成一个mn阶上三角阵；triu(A,k)将从矩阵A中提取主对角线第|k|条对角线之上的上三角部分构成一个mn阶上三角阵。注意：这里的k与diag(A,k)的用法类似，当k＞0，则该对角线位于主对角线的上方第k条；当k＜0，该对角线位于主对角线的下方第|k|条；当k=0，则等同于triu(A)【例4】试分别用triu(A)、triu(A,1)和、triu(A,-1)从矩阵A提取相应的上三角部分构成上三角阵B、C和D。MATLAB程序如下：A=[123;456;789;987];%一个已知的43阶矩阵A%构成各种情况的上三角阵B、C和DB=triu(A)B=123056009000C=triu(A,1)D=triu(A,-1)9.下三角阵：使用格式为tril(A)、tril(A,k)tril的功能是从矩阵A中提取下三角部分构成下三角阵。用法与triu相同。10.空矩阵在MATLAB里，把行数、列数为零的矩阵定义为空矩阵。空矩阵在数学意义上讲是空的，但在MATLAB里确是很有用的。例如A=[0.10.20.3;0.40.50.6];B=find(A1.0)B=[]这里[]是空矩阵的符号，B=find(A1.0)表示列出矩阵A中值大于1.0的元素的序号。当不能满足括号中的条件时，返回空矩阵。另外，也可以将空矩阵赋给一个变量，如：B=[]B=[]Hilbert矩阵及逆Hilbert矩阵invhilb(n)专门求希尔伯特矩阵的逆的函数Vandermonde矩阵【例4-4】魔方矩阵魔方矩阵：每行、每列及两条对角线上的元素和都相等。对于n阶魔方阵，其元素由1,2,3,…,n2共n2个整数组成。函数magic(n)生成一个n阶魔方阵。托普利兹矩阵托普利兹(Toeplitz)矩阵除第一行第一列外，其他每个元素都与左上角的元素相同。1678216732164321toeplitz(x,y)：生成托普利兹矩阵的函数，它生成一个以x为第一列，y为第一行的托普利兹矩阵。这里x,y均为向量，两者不必等长。toeplitz(x)：用向量x生成一个对称的托普利兹矩阵。例如T=toeplitz(1:6)toeplitz([1234],[1678])伴随矩阵p(x)=0的(所有)根是A的(所有)特征值.11011101.()...)..(nnnnnnnnnnnnaaaxpxaxaxaxxxqxaaaaaa（）1231000000000000000011100nnnnnnnnaaaaaaaaaaA2(1)()nqAEp(x)是A的特征多项式110[,,...,,];compan();nnPaaaaP例如，求多项式的x3-7x+6的伴随矩阵，可使用命令：p=[1,0,-7,6];compan(p)帕斯卡矩阵例3.5求(x+y)5的展开式。命令：pascal(6)矩阵次对角线上的元素1,5,10,10,5,1即为展开式的系数。111111123456136101521141020355615153570126162156126252二次项(x+y)n展开后的系数随n的增大组成一个三角形数表，称为杨辉三角形。由杨辉三角形生成的矩阵称为帕斯卡(Pascal)矩阵。函数pascal(n)生成一个n阶帕斯卡矩阵。二、矩阵的特征值与特征向量对于N×N阶方阵A，所谓A的特征值问题是：求数λ和N维非零向量x（通常为复数），使之满足下式：Ax=λx则称λ为矩阵A的一个特征值（特征根），而非零向量x为矩阵A的特征值λ所对应的特征向量。对一般的N×N阶方阵A，其特征值通常为复数，若A为实对称矩阵，则A的特征值为实数。MATLAB提供的内部函数eig可以用来计算特征值与特征向量。eig函数的使用格式有五种，其中常见的有E=eig(A)、[V,D]=eig(A)和[V,D]=eig(A,’nobalance’)三种，另外两种格式用来计算矩阵的广义特征值与特征向量：E=eig(A,B)和[V,D]=eig(A,B)。(1)E=eig(A)：由eig(A)返回方阵A的N个特征值，构成向量E；(2)[V,D]=eig(A)：由eig(A)返回方阵A的N个特征值，构成NN阶对角阵D，其对角线上的N个元素即为相应的特征值，同时将返回相应的特征向量赋予NxN阶方阵V的对应列，且A、V、D满足AV=VD；(3)[V,D]=eig(A,’nobalance’)：本格式的功能与格式(2)一样，只是格式(2)是先对A作相似变换(balance)，然后再求其特征值与相应的特征向量；而本格式则事先不作相似变换；(4)E=eig(A,B)：由eig(A,B)返回N×N阶方阵A和B的N个广义特征值，构成向量E。(5)[V,D]=eig(A,B)：由eig(A,B)返回方阵A和B的N个广义特征值，构成NN阶对角阵D，其对角线上的N个元素即为相应的广义特征值，同时将返回相应的特征向量构成N×N阶满秩矩阵，且满足AV=BVD。【例5】试用格式(1)求下列对称矩阵A的特征值；用格式(2)求A的特征值和相应的特征向量，且验证之。A=[1.00001.00000.50001.00001.00000.25000.50000.25002.0000];执行eig(A)将直接获得对称矩阵A的三个实特征值：eig(A)ans=-0.01661.48012.5365而下列命令则将其三个实特征值作为向量赋予变量E：E=eig(A)E=-0.01661.48012.5365矩阵的秩与迹1．矩阵的秩矩阵线性无关的行数或列数称为矩阵的秩。求矩阵秩的函数：rank(A)。2．矩阵的迹矩阵的迹等于矩阵的对角线元素之和，也等于矩阵的特征值之和。（证明：《高等代数》或《矩阵论》）求矩阵的迹的函数是trace(A)。求向量范数的函数：(1)norm(V)或norm(V,2)：计算向量V的2—范数。(2)norm(V,1)：计算向量V的1—范数。(3)norm(V,inf)：计算向量V的∞—范数。(1)2-范数(2)1-范数(3)∞-范数221Vniiv1maxiinVv1．向量的3种常用范数及其计算函数11Vniiv123(,,,,)nVvvvv三、行列式的值MATLAB提供的内部函数det用来计算矩阵的行列式的值。设矩阵A为一方阵(必须是方阵)，求矩阵A的行列式值的格式为：det(A)。注意：本函数同样能计算通过构造出的稀疏矩阵的行列式的值。关于如何构造稀疏矩阵，将在本章最后一节介绍。【例6】利用随机函数产生一个三阶方阵A，然后计算方阵之行列式的值。A=rand(3)A=0.95010.48600.45650.23110.89130.01850.60680.76210.8214det(A)ans=0.4289四、矩阵求逆及其线性代数方程组求解若方阵A,B满足等式A*B=B*A=I(I为单位矩阵)则称A为B的逆矩阵，或称B为A的逆矩阵。这时A，B都称为可逆矩阵(或非奇异矩阵、或满秩矩阵)，否则称为不可逆矩阵(或奇异矩阵、或降秩矩阵)。矩阵求逆【例7】试用inv函数求方阵A的逆阵A-1赋值给B，且验证A与A-1是互逆的。A=[1-11;5-43;211];B=inv(A)B=-1.40000.40000.20000.2000-0.20000.40002.6000-0.60000.2000A*Bans=1.00000.0