Matlab第4章

第4章线性代数中的数值计算【本章学习目标】●掌握生成特殊矩阵的方法。●掌握矩阵分析的方法。●掌握求解线性方程组的各种方法。●了解矩阵的稀疏存储方式。4.1特殊矩阵的生成4.1.1通用的特殊矩阵●zeros函数：产生全0矩阵，即零矩阵。●ones函数：产生全1矩阵，即幺矩阵。●eye函数：产生单位矩阵，即对角线上的元素为1、其余元素为0的矩阵。●rand函数：产生0～1均匀分布的随机矩阵。●randn函数：产生均值为0、方差为1的标准正态分布随机矩阵。这几个函数的调用格式相似，如果这个函数的参数只是一个，那么MATLAB将会创建一个方阵，行数和列数均为这个参数；如果这个函数的参数有两个，那么第一个参数代表行数，第二个参数代表列数。下面以产生零矩阵的zeros函数为例进行说明。zeros函数的调用格式如下。●zeros(m)：产生m×m零矩阵。●zeros(m,n)：产生m×n零矩阵。当m=n时，等同于zeros(m)。●zeros(size(A))：产生与矩阵A同样大小的零矩阵。【例4.1】分别建立3×3、3×2和与矩阵A同样大小的零矩阵。（1）建立一个3×3的零矩阵。zeros(3)ans=000000000（2）建立一个2×3的零矩阵。zeros(2,3)（3）设A为2×3矩阵，则可以用zeros(size(A))建立一个与矩阵A同样大小的零矩阵。A=[123;456];%产生一个2×3阶矩阵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+(b−a)xi计算即可。产生均值为0、方差为1的标准正态分布随机矩阵使用randn函数，假设已经得到了一组标准正态分布随机数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】将101～125等25个数填入一个5行5列的表格中，使其每行、每列及对角线的和均为565。一个5阶魔方矩阵的每行、每列及对角线的和均为65，对其每个元素都加100后这些和变为565。完成其功能的命令如下：M=100+magic(5)M=1171241011081151231051071141161041061131201221101121191211031111181251021092．范得蒙矩阵范得蒙（Vandermonde）矩阵的最后一列全为1，倒数第二列为一个指定的向量，其他各列是其后列与倒数第二列的点乘积。可以用一个指定向量生成一个范得蒙矩阵。在MATLAB中，函数vander(V)生成以向量V为基础向量的范得蒙矩阵。3．希尔伯特矩阵希尔伯特（Hilbert）矩阵是一种数学变换矩阵，它的每个元素hij=1/(i+j−1)。在MATLAB中，生成希尔伯特矩阵的函数是hilb(n)。专门求希尔伯特矩阵的逆的函数invhilb(n)4．托普利兹矩阵托普利兹（Toeplitz）矩阵除第一行第一列外，其他每个元素都与左上角的元素相同。生成托普利兹矩阵的函数是toeplitz(x,y)，它生成一个以x为第一列、y为第一行的托普利兹矩阵。这里x、y均为向量，两者不必等长。toeplitz(x)用向量x生成一个对称的托普利兹矩阵。5．伴随矩阵生成伴随矩阵的函数是compan(p)，其中p是一个多项式的系数向量，高次幂系数排在前，低次幂排在后。例如，为了求多项式的x3−7x+6的伴随矩阵，可使用如下命令：p=[1,0,-7,6];compan(p)ans=07-61000106．帕斯卡矩阵我们知道，二次项(x+y)n展开后的系数随n的增大组成一个三角形表，称为杨辉三角形。由杨辉三角形表组成的矩阵称为帕斯卡（Pascal）矩阵，它的元素p1j=1，pi1=1，pij=pi−1,j−1+pi−1,j（i1，j1）。函数pascal(n)生成一个n阶帕斯卡矩阵。【例4.5】求(x+y)4的展开式。在MATLAB命令窗口，输入命令：pascal(5)ans=111111234513610151410203515153570矩阵次对角线上的元素1，4，6，4，1即为展开式的系数，即(x+y)4=x4+4x3y+6x2y2+4xy3+y44.2.1矩阵结构变换1．对角阵只有对角线上有非0元素的矩阵称为对角矩阵，对角线上的元素都为1的对角矩阵称为单位矩阵。（1）提取矩阵的对角线元素设A为m×n矩阵，函数diag(A)用于提取矩阵A主对角线元素，产生一个具有min(m,n)个元素的列向量。例如：A=[1,2,3;4,5,6]A=123456D=diag(A)D=15diag(A,k)提取第k条对角线的元素。主对角线为第0条对角线；与主对角线平行，往上为第1条，第2条，…，第n条对角线，往下为第−1条，第−2条，…，第−n条对角线。（2）构造对角矩阵设V为具有m个元素的向量，diag(V,k)的功能是产生一个n×n（n=m+|k|）对角阵，其第k条对角线的元素即为向量V的元素。例如：diag(1:3,-1)ans=0000100002000030省略k时，相当于k为0，其主对角线元素即为向量V的元素。【例4.6】先建立5×5矩阵A，然后将A的第一行元素乘以1，第二行乘以2，…，第五行乘以5。用一个对角矩阵左乘一个矩阵时，相当于用对角阵的第一个元素乘以该矩阵的第一行，用对角阵的第二个元素乘以该矩阵的第二行……依此类推，因此，只需按要求构造一个对角矩阵D，并用D左乘A即可。命令如下：A=[1:5;2:6;3:7;4:8;5:9]D=diag(1:5);D*A%用D左乘A,对A的每行乘以一个指定常数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列的矩阵经过转置运算后，变成一个n行m列的矩阵。MATLAB中，转置运算符是单撇号(')。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=38853731849-20rot90(A,4)ans=93738-231808455．矩阵的翻转矩阵的翻转分左右翻转和上下翻转。对矩阵实施左右翻转是将原矩阵的第一列和最后一列调换，第二列和倒数第二列调换……依此类推。对矩阵A实施左右翻转的函数是fliplr(A)。对矩阵A实施上下翻转的函数是flipud(A)。4.2.2矩阵求值1．方阵的行列式值把一个方阵看作一个行列式，并对其按行列式的规则求值，这个值就称为矩阵所对应的行列式的值。在MATLAB中，求方阵A所对应的行列式的值的函数是det(A)。例如：A=[1:3;2:-1:0;12,5,9]A=1232101259B=det(A)B=-332．矩阵的秩与迹（1）矩阵的秩rank(A)（2）矩阵的迹矩阵的迹即矩阵的对角线元素之和。trace(A)。3．向量和矩阵的范数设向量V=(v1,v2,…,vn)，向量的3种常用范数如下。1-范数：11niiVv。2-范数：221niiVv。?-范数：1max{}iinVv≤≤。设A是一个m×n的矩阵，矩阵的3种常用范数如下。1-范数：111max{}mijjniAa≤≤。2-范数：12A，其中λ1为A'A最大特征值。∞-范数：11max{}nijimjAa≤≤。函数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−范数下的条件数。●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的特征值和特征向量，而格式3直接求矩阵A的特征值和特征向量。一个矩阵的特征向量有无穷多个，eig函数只找出其中的n个，A的其他特征向量，均可由这n个特征向量的线性组合表示。【例4.7】用求特征值的方法解方程：3x5−7x4+5x2+2x−18=0先构造与方程对应的多项式的伴随矩阵A，再求A的特征值。A的特征值即为方程的根。命令如下：p=[3,-7,0,5,2,-18];A=compan(p);%构造伴随矩阵Ax1=eig(A)%求A的特征值x1=2.18371.0000+1.0000i1.0000-1.0000i-0.9252+0.7197i-0.9252-0.7197ix2=roots(p)%直接求多项式p的零点x2=2.18371.0000+1.0000i1.0000-1.0000i-0.9252+0.7197i-0.9252-0.7197i4.3线性方程组求解4.3.1矩阵求逆及线性代数方程组求解1．矩阵求逆inv(A)函数用于计算方阵的逆矩阵。2．利用矩阵求逆方法解线性方程组【例4.9】利用矩阵求逆方法解线性方程组A=[1,-2,3;3,-1,5;2,1,5];b=[1;2;3];x=inv(A)*bx=-0.33330.33330.6667123123123231352253xxxxxxxxx4.3.2利用左除运算符求解线性方程组对于线性方程组Ax=b，可以利用左除运算符“\”求解：x=A\b【例4.10】用左除运算符求解下列相同系数矩阵的两个线性代数方程组的解。（1）123111254332111xxx（2