第三章_MATLAB矩阵分析与处理

1、 第3章 MATLAB矩阵分析与处理 特殊矩阵 矩阵结构变换 矩阵求逆与线性方程组求解 矩阵求值 3.1 特殊矩阵 3.1.1 通用的特殊矩阵 常用的产生通用特殊矩阵的函数有： zeros：产生全0矩阵(零矩阵)。 ones：产生全1矩阵(幺矩阵)。 eye：产生单位矩阵。 rand：产生01间均匀分布的随机矩阵。 randn：产生均值为0，方差为1的标准正态分布随机 矩阵。 这几个函数的调用格式相似，下面以产生零矩 阵的zeros函数为例进行说明。其调用格式为： zeros(m): 产生mm零矩阵 zeros(m,n) 产生mn零矩阵 zeros(size(A) 产生与矩阵A同样大小的零矩阵

2、。 例3.1 分别建立33、32和与矩阵A同样大小的零矩阵 (1) 建立一个33零矩阵。 zeros(3) ans= 0 0 0 0 0 0 0 0 0 (2) 建立一个32零矩阵。 zeros(3,2) ans= 0 0 0 0 0 0 (3) 设A为23矩阵，则可以用zeros(size(A)建立一个与 矩阵A同样大小零矩阵。 A=1 2 3;4 5 6; %产生一个23阶矩阵A zeros(size(A) %产生一个与矩阵A同样大小的零矩阵 ans= 0 0 0 0 0 0 例3.2 建立随机矩阵： (1) 在区间20,50内均匀分布的5阶随机矩阵。 (2) 均值为0.6、方差为0.1的

3、5阶正态分布随机矩阵。 rand：产生01间均匀分布的随机矩阵。 要得到a,b区间上均匀分布的随机数，需用yi=a+(b-a)xi randn：产生均值为0，方差为1的标准正态分布随机矩阵。 命令如下： x=20+(50-20)*rand(5) x = 48.5039 42.8629 38.4630 32.1712 21.7367 26.9342 33.6940 43.7581 48.0641 30.5860 38.2053 20.5551 47.6544 47.5071 44.3950 34.5795 44.6422 42.1462 32.3081 20.2958 46.7390 33.34

4、11 25.2880 46.8095 24.1667 y=0.6+sqrt(0.1)*randn(5) y = 0.8713 0.4735 0.8114 0.0927 0.7672 0.9966 0.8182 0.9766 0.6814 0.6694 0.0960 0.8579 0.2197 0.2659 0.3085 0.1443 0.8251 0.5937 1.0475 -0.0864 0.7806 1.0080 0.5504 0.3454 0.5813 3.1.2 用于专门学科的特殊矩阵 (1) 魔方矩阵 魔方矩阵有一个有趣的性质，其每行、每 列及两条对角线上的元素和都相等。对于n 阶魔

5、方阵，其元素由1,2,3,n2共n2个整数 组成。MATLAB提供了求魔方矩阵的函数 magic(n)，其功能是生成一个n阶魔方阵。 magic(3) ans = 8 1 6 3 5 7 4 9 2 例3.3 将101-125等25个数填入一个5行5列的表格 中，使其每行每列及对角线的和均为565。 一个5姐魔方矩阵的每行、每列及对角线的和 均为65，对其每个元素都加100后，这些和变成565. magic(5)magic(5) ans =ans = 17 24 1 8 15 17 24 1 8 15 23 5 7 14 16 23 5 7 14 16 4 6 13 20 22 4 6 13

6、20 22 10 12 19 21 3 10 12 19 21 3 11 18 25 2 9 11 18 25 2 9 M=100+magic(5)M=100+magic(5) M =M = 117 124 101 108 115 117 124 101 108 115 123 105 107 114 116 123 105 107 114 116 104 106 113 120 122 104 106 113 120 122 110 112 119 121 103 110 112 119 121 103 111 118 125 102 109 111 118 125 102 109 (2)

7、范得蒙德矩阵 范得蒙德(Vandermonde)矩阵最后一列全为1， 倒数第二列为一个指定的向量，其他各列是其后列 与倒数第二列的点乘积。可以用一个指定向量生成 一个范得蒙矩阵。在MATLAB中，函数vander(V)生 成以向量V为基础向量的范得蒙矩阵。例如， A=vander(1;2;3;5)即可得到上述范得蒙矩阵。 A =A = 1 1 1 1 1 1 1 1 8 4 2 1 8 4 2 1 27 9 3 1 27 9 3 1 125 25 5 1 125 25 5 1 在MATLAB中，生成希尔伯特矩阵的函数是hilb(n). 使用一般方法求逆会因为原始数据的微小扰动而产 生不可靠的计

8、算结果。MATLAB中，有一个专门求希 尔伯特矩阵的逆的函数invhilb(n)，其功能是求n阶 的希尔伯特矩阵的逆矩阵。 (3) (3) 希尔伯特矩阵希尔伯特矩阵 1 (). 1 ij Hilberth ij 希尔伯特矩阵的每个元素 例3.4 求4阶希尔伯特矩阵及其逆矩阵。 命令如下： hilb(4)hilb(4) ans =ans = 1 1/2 1/3 1/4 1 1/2 1/3 1/4 1/2 1/3 1/4 1/5 1/2 1/3 1/4 1/5 1/3 1/4 1/5 1/6 1/3 1/4 1/5 1/6 1/4 1/5 1/6 1/7 1/4 1/5 1/6 1/7 invhi

9、lb(4)invhilb(4) ans =ans = 16 -120 240 -140 16 -120 240 -140 -120 1200 -2700 1680 -120 1200 -2700 1680 240 -2700 6480 -4200 240 -2700 6480 -4200 -140 1680 -4200 2800 -140 1680 -4200 2800 (4) 托普利兹矩阵 托普利兹(Toeplitz)矩阵除第一行第一列外，其他每个元素 都与左上角的元素相同。生成托普利兹矩阵的函数是 toeplitz(x,y)，它生成一个以x为第一列，y为第一行的托普利 兹矩阵。这里x, y

10、均为向量，两者不必等长,toeplitz(x)用向量x 生成一个对称的托普利兹矩阵。如:T=toeplitz(1:6) T=toeplitz(1:4) T = 1 2 3 4 2 1 2 3 3 2 1 2 4 3 2 1 (5) 伴随矩阵 1 110 12301 ( )( ) 10000 01000 000100 000010 ( )( ) ( ) nn nn nnn nnnnn p xp xa xaxa xa aaaaa aaaaa A p xp xA p x 设多项式为： 称矩阵： 为多项式的伴随矩阵。成为 的特征多项式，方程 的根成为A的特征值。 MATLAB生成伴随矩阵的函数是com

11、pan(p)，其 中p是一个多项式的系数向量，高次幂系数排在前， 低次幂排在后。例如，为了求多项式的x3-7x+6的伴 随矩阵，可使用命令： p=1,0,-7,6; p=1,0,-7,6; compan(p) compan(p) ans =ans = 0 7 -6 0 7 -6 1 0 0 1 0 0 0 1 0 0 1 0 (6) 帕斯卡矩阵 我们知道，二次项(x+y)n展开后的系数随n 的增大组成一个三角形表，称为杨辉三角 形。由杨辉三角形表组成的矩阵称为帕斯 卡(Pascal)矩阵。函数pascal(n)生成一个n阶 帕斯卡矩阵。 例3.5 求(x+y)5的展开式。 在MATLAB命令窗

12、口，输入命令： pascal(6) 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 2 矩阵次对角线上的元素1,5,10,10,5,1即为展开式 的系数。 3.2 矩阵结构变换 3.2.1 对角阵与三角阵 1对角阵 只有对角线上有非0元素的矩阵称为对角矩阵， 对角线上的元素相等的对角矩阵称为数量矩阵，对 角线上的元素都为1的对角矩阵称为单位矩阵。 矩阵对角线有很多性质，如转置矩阵时对角线元 素不变，相似变换时对角线的和（称为矩阵的迹） 不变

13、等。 (1) 提取矩阵的对角线元素 设A为mn矩阵，diag(A)函数用于提取矩阵A主对角线 元素，产生一个具有min(m,n)个元素的列向量。 A=1 2 3;4 5 6; D=diag(A) D = 1 5 diag(A)函数还有一种形式diag(A,k)，其功能是提取第 k条对角线的元素。与主对角线平行，往上为第1条，第2条， ，第n条对角线，往下为第-1条，第-2条，第-n条对 角线。主对角线为第0条对角线。例如对上面建立的A矩阵， 提取主对角线两侧对角线的元素，命令如下： D1=diag(A,1)D1=diag(A,1) D1 =D1 = 2 2 6 6 D2=diag(A,-1)D

14、2=diag(A,-1) D2 =D2 = 4 4 (2) (2) 构造对角矩阵构造对角矩阵 设设V V为具有为具有m m个元素的向量，个元素的向量，diag(V)diag(V)将产生一将产生一 个个m mm m对角矩阵，其主对角线元素即为向量对角矩阵，其主对角线元素即为向量V V的元素的元素 diag(1,2,-1,4)diag(1,2,-1,4) ans =ans = 1 0 0 0 1 0 0 0 0 2 0 0 0 2 0 0 0 0 -1 0 0 0 -1 0 0 0 0 4 0 0 0 4 diag(V)函数也有另一种形式diag(V,k)，其功 能是产生一个nn(n=m+|k|)

15、对角阵，其第k 条对角线的元素即为向量V的元素。 diag(1:3,-1) ans = 0 0 0 0 1 0 0 0 0 2 0 0 0 0 3 0 例3.6 先建立55矩阵A，然后将A的第一行元 素乘以1，第二行乘以2，第五行乘以5。 A=17,0,1,0,15;23,5,7,14,16;4,0,13,0,22;10,12,19,21,3;. 1,18,25,2,19; D=diag(1:5); D*A %用D左乘A，对A的每行乘以一个指定常数 ans = 17 0 1 0 15 46 10 14 28 32 12 0 39 0 66 40 48 76 84 12 55 90 125 10

16、 95 对对A A的每列元素乘的每列元素乘 以同一个数，可以同一个数，可 以用一个对角阵以用一个对角阵 右乘右乘A.A. 2三角阵 三角阵又进一步分为上三角阵和下三角阵，所 谓上三角阵，即矩阵的对角线以下的元素全为0的一 种矩阵，而下三角阵则是对角线以上的元素全为0的 一种矩阵。 (1) 上三角矩阵 与矩阵A对应的上三角阵B是与A同型的一个矩阵，并且B 的对角线以上（含对角线）和A对应相等，而对角线以下的 元素等于0。求矩阵A的上三角阵的MATLAB函数是triu(A)。 例如，提取矩阵A的上三角元素，形成新的矩阵B。 A=7,13,-28;2,-9,8;0,34,5; B=triu(A) B

17、 = 7 13 -28 0 -9 8 0 0 5 triu(A)函数也有另一种形式triu(A,k)，其功能是求矩 阵A的第k条对角线以上的元素。例如，提取矩阵A的第2条对 角线以上的元素，形成新的矩阵B。 A=1,32,1,0,5;3,5,17,4,16;4,0,13,0,42;70,11,9,21,3;11,63,5,2,99; B=triu(A,2) B = 0 0 1 0 5 0 0 0 4 16 0 0 0 0 42 0 0 0 0 0 0 0 0 0 0 (2) 下三角矩阵 在MATLAB中，提取矩阵A的下三角矩阵的 函数是 tril(A)和tril(A,k)，其用法与提取上三角矩

18、阵的 函数triu(A)和triu(A,k)完全相同。 3.2.2 矩阵的转置与旋转 1矩阵的转置 所谓转置，即把源矩阵的第一行变成目标矩阵 的第一列，第二行变成第二列，依此类推。 显 然，一个m行n列的矩阵经过转置运算后，变成 一个 n行m列的矩阵。转置运算符是单撇号()。 A=71,3,-8;2,-9,8;0,4,5; B=A B = 71 2 0 3 -9 4 -8 8 5 2矩阵的旋转 在MATLAB中，可以很方便地以90。为单位对矩阵 按逆时针方向旋转：。利用函数rot90(A,k)将矩阵A旋 转90的k倍，当k为1时可省略。例如，将A按逆时针 旋转90。，命令如下 A=57,19,

19、38;-2,31,8;0,84,5;A=57,19,38;-2,31,8;0,84,5; B=rot90(A)B=rot90(A) B =B = 38 8 5 38 8 5 19 31 84 19 31 84 57 -2 0 57 -2 0 rot90(A,4)rot90(A,4) ans =ans = 57 19 38 57 19 38 -2 31 8 -2 31 8 0 84 5 0 84 5 3矩阵的左右翻转 对矩阵实施左右翻转是将原矩阵的第一列和 最后一列调换，第二列和倒数第二列调 换，依次类推。MATLAB对矩阵A实施左 右翻转的函数是fliplr(A) A=14,-9,8;-2,8

20、1,8;-2,4,0 A = 14 -9 8 -2 81 8 -2 4 0 B=fliplr(A) B = 8 -9 14 8 81 -2 0 4 -2 4矩阵的上下翻转 与矩阵的左右翻转类似，矩阵的上下翻转 是将原矩阵的第一行与最后一行调换，第 二行与倒数第二行调换，依次类推。 MATLAB对矩阵A实施上下翻转的函数是 flipud(A)。 3.3 矩阵求逆与线性方程组求解 3.3.1 矩阵的逆与伪逆 对于一个方阵A，如果存在一个与其同阶的方 阵B，使得： AB=BA=I (I为单位矩阵) 则称B为A的逆矩阵，当然，A也是B的逆矩阵。 求一个矩阵的逆是一件非常烦琐的工作，容易 出错，但在MA

21、TLAB中，求一个矩阵的逆非常容易。 求方阵A的逆矩阵可调用函数inv(A)。 例3.7 求方阵A的逆矩阵，且验证A与A-1是否是互逆 的。 A=1,-1,1;5,-4,3;2,1,1; B=inv(A); A*B ans = 1.0000 0.0000 0.0000 -0.0000 1.0000 0.0000 -0.0000 0.0000 1.0000 B*A ans = 1.0000 0.0000 -0.0000 -0.0000 1.0000 0.0000 0.0000 -0.0000 1.0000 上述计算中可上述计算中可 见：见：AB=BAAB=BA 即：即：AAAA-1 -1=A =

22、A-1 -1A, A, 故故A A与与A-1A-1是互是互 逆的。逆的。 如果矩阵A不是一个方阵，或者A是一个非满秩 的方阵时，矩阵A没有逆矩阵，但可以找到一 个与A的转置矩阵A同型的矩阵B，使得： ABA=A BAB=B 此时称矩阵B为矩阵A的伪逆，也称为广义逆 矩阵。在MATLAB中，求一个矩阵伪逆的函数 是： pinv(A) A=3,1,1,1;1,3,1,1;1,1,3,1; B=pinv(A) B = 0.3929 -0.1071 -0.1071 -0.1071 0.3929 -0.1071 -0.1071 -0.1071 0.3929 0.0357 0.0357 0.0357 若A

23、是一个奇异矩阵，无一般意义上的逆矩阵， 但可以求A的伪逆矩阵。例如： A=0,0,0;0,1,0;0,0,1; pinv(A) ans = 0 0 0 0 1 0 0 0 1 本例中，A的伪逆矩阵和A相等，这是一个巧 合。一般说来，矩阵的伪逆矩阵和自身是不同的。 将包含n个未知数，由n个方程构成的线性方程 组表示成： 3.2.2 用矩阵求逆方法求解线性方程组 11112211 21122221 11221 nn nn nnnnn a xa xa xb a xa xa xb a xa xa xb 1112111 2122222 12 , n n nnnnnn aaaxb aaaxb Axb a

24、Axb aaxb 其矩阵表示为： 其： 中 在线性方程组Ax=b两边各左乘A-1,有： A-1 Ax= A-1 b 由于A-1 A=I，故得： x= A-1 b 所以，利用系数矩阵A的逆矩阵，可以求解线 性方程组。 3.8 235 492 8276 xyz xyz xyz 例用求逆矩阵的方法解线性方程组。 命令如下：命令如下： A=1,2,3;1,4,9;1,8,27;A=1,2,3;1,4,9;1,8,27; b=5,-2,6b=5,-2,6; ; x=inv(A)x=inv(A)* *b b x =x = 23.0000 23.0000 -14.5000 -14.5000 3.6667 3

25、.6667 也可以运用左除运算也可以运用左除运算 符符“ ”求解线性代数方程求解线性代数方程 组。组。 A=1,2,3;1,4,9;1,8,27;A=1,2,3;1,4,9;1,8,27; b=5,-2,6b=5,-2,6; ; x=Abx=Ab 3.4 矩阵求值 3.4.1方阵的行列式 把一个方阵看作一个行列式，并对其按行列式 的规则求值，这个值就称为矩阵所对应的行列式的 值。在MATLAB中，求方阵A所对应的行列式的值的函 数是det(A)。 A=rand(5) A = 0.9501 0.7621 0.6154 0.4057 0.0579 0.2311 0.4565 0.7919 0.93

26、55 0.3529 0.6068 0.0185 0.9218 0.9169 0.8132 0.4860 0.8214 0.7382 0.4103 0.0099 0.8913 0.4447 0.1763 0.8936 0.1389 B=det(A) B = -0.0071 1矩阵的秩 矩阵线性无关的行数与列数称为矩阵的秩。在 MATLAB中，求矩阵秩的函数是rank(A)。 3.4.2 矩阵的秩与迹 1 122 2 1 , (1,2, 0 ) pp p i k xk x xx k p xk x ip 对于一组向量若存在一组不全为零 ，使： 成立，则称 个向量线性相关，否则称线性相关。 A=2,2

27、,-1,1;4,3,-1,2;8,5,-3,5;3,3,-2,2; r=rank(A) r = 4 这说明A是一个满秩矩阵。 2矩阵的迹 矩阵的迹等于矩阵的对角线元素之和，也等 于矩阵的特征值之和。在MATLAB中，求矩阵 的迹的函数是trace(A)。 A=2,2,3;4,5,-6;7,8,9; trace(A) ans = 16 矩阵或向量的范数用来度量矩阵或向量在某种 意义下的长度。范数有多种方法定义，其定义不 同，范数值也就不同。 3.4.3 向量和矩阵的范数 12 2 2 1 1 1 1 ,3 (1)2 (2)1 (3)m 1 x 3 a n n i i n i i i i n Vv

28、 vv Vv Vv Vv 设向量下面讨论向量的 种范数： 范数 范数 、向量的 种常见范数及其 范数 计 算函数 在MATLAB中，求这3种向量范数的函数分别为： (1) norm(V)或norm(V,2)：计算向量V的2-范数 (2) norm(V,1)：计算向量V的1-范数。 (3) norm(V,inf)：计算向量V的-范数。 例如：V=-1,1/2,1; V1=norm(V,1) %求V的1-范数 V1 = 2.5000 V2=norm(V) %求V的2-范数 V2 = 1.5000 Vinf=norm(V,inf) %求V的-范数 Vinf = 1 2矩阵的范数及其计算函数 设A是一

29、个mn的矩阵，V是一个含有n个元素的 列向量，定义： 因为A是一个mn的矩阵，而V是一个含有n个元 素的列向量。在前面已经定义了3种不同的向量范 数，按照上式也可以定义3种矩阵范数，这样定义的 矩阵范数A称为A从属于向量的范数。 max|,|,| | 1|A|=A Vv MATLAB提供了求3种矩阵范数的函数，其函数调 用格式与求向量的范数的函数完全相同。 1 11 1 11 2 1 2 1 1 1 1 3( ) maxmax() max maxmax() ij m ij Vj n i V n ij Vi n j a A A Va A VA A A V A Aa A 从属于 种向量范数的矩阵范

30、数计算公式是是 矩阵 的元素 ： 按列加 ，其中 为最大特征值 按行加 A=17,0,1,0,15;23,5,7,14,16;4,0,13,0,22;10,12,19,2 1,3;11,18,25,2,19; a1=norm(A,1) %求A的1-范数 a1 = 75 a2=norm(A,2) %求A的2-范数 a2 = 59.3617 ainf=norm(A,inf) %求A的-范数 ainf = 75 3.4.4 矩阵的条件数 在求解线性方程组Ax=b时，一般认为：系数矩 阵A中个别元素的微小扰动不会引起解向量的很大变 化。这样的假设在工程应用中非常重要，因为一般 系数矩阵是由实验数据获得

31、的，并非精确解，但与 精确解误差不大。基于上述假设可以得到如下结 论：当参与运算的系数与实际精确解误差很小时， 所获得的解与问题的精确解误差也很小。 对于有的系数矩阵，个别元素的微小扰动会引 起解的很大变化，在计算数学中，称这种矩阵是病 态矩阵。 而称解不因系数矩阵的微小扰动而发生的大的 变化的矩阵为良性矩阵。 当然良性与病态也是相对的，需要一个参数来描 述，条件数就是用来描述矩阵的这种性能的一个参 数。 矩阵A的条件数等于A的范数与A的逆矩阵的范数 的乘积，即cond(A)=AA-1。这样定义的条 件数总是大于1的。条件数越接近于1，矩阵的性能 越好，反之，矩阵的性能越差。 在MATLAB中

32、，计算矩阵A的3种条件数的函数是： (1) cond(A,1) 计算A的1-范数下的条件数。 即：cond(A,1)=A1A-11 (2) cond(A)或cond(A,2) 计算A的2-范数数下 的条件数。 即： cond(A)=A2A-12 (3) cond(A,inf) 计算A的-范数下的条件数。 即：cond(A,inf)=cond(A)=AA-1 例如 : A=2,2,3;4,5,-6;7,8,9; C1=cond(A) C1 = 87.9754 B=2,-5,4;1,5,-2;-1,2,4; C2=cond(B) C2 = 3.7515 矩阵B的条件数比矩阵A的条件数更接近于1，因

33、此，矩 阵B的性能要好于矩阵A。 3.5 矩阵的特征值与特征向量 nAA A A 对于 阶方阵 ，求数 和向量 ，使得等式 成立。满足等式的数 称为矩阵 的，而特向量 称 为 的 征值 特征向量。 0 0 0 AA AI I AI 方程和是两个等价方程，要 使方程有非零解 ，必须使其系数行列式 为零，即。 0 AIn AInn An 行列式是一个关于 的阶多项式，因而 方程是一个次方程，有 个根（含重根）， 就是矩阵 的 个特征值，每一个特征值对应无穷多个 特征向量。矩阵的特征值问题有确定解，但特征向量 没有确定解。 在在MATLABMATLAB中，计算矩阵中，计算矩阵A A的特征值和特征向量

34、的的特征值和特征向量的 函数是函数是eig(A)eig(A)，常用的调用格式有，常用的调用格式有3 3种：种： (1)E=eig(A)(1)E=eig(A)：求矩阵求矩阵A A的全部特征值，构成向量的全部特征值，构成向量E E。 (2) V,D=eig(A)：求矩阵A的全部特征值，构成 对角阵D，并求A的特征向量构成V的列向量。 (3) V,D=eig(A,nobalance):与第二种格式类 似，但第二种格式中先对A做相似变换后，求矩阵A 的特征值和特征向量，而格式3直接求矩阵A的特征 值和特征向量。 一个矩阵的特征向量有无穷多个，eig函数只找 出其中的n个，A的其他特征向量均可由n个特征

35、向量 的线性组合表示。 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.7103 D = -0.0166 0 0 0 1.4801 0 0 0 2.5365 求得的3个特征值是-0.0166，1.4801，2.5365 各特征值对应的特征向量为V的各列的向量。验 证结 果，AV和VD的值均为： -0.0120 0.6576 1.3481 0.0114 0.8320 1.1705 0.0016 -1.0325 1.8018 例3.9

36、 用求特征值的方法解方程。 3x5-7x4+5x2+2x-18=0 先构造与方程对应的多项式的伴随矩阵A，再 求A的 特征值。A的特征值即为方程的根。命令如下： p=3,-7,0,5,2,-18; A=compan(p); %A的伴随矩阵 x1=eig(A) %求A的特征值 x1 = 2.1837 1.0000 + 1.0000i 1.0000 - 1.0000i -0.9252 + 0.7197i -0.9252 - 0.7197i x2=roots(p) %直接求多项式p的零点 x2 = 2.1837 1.0000 + 1.0000i 1.0000 - 1.0000i -0.9252 + 0.7197i -0.9252 - 0.7197i 可以看出，两种方法求得的方程的