线性代数中的数值计算.ppt

第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的标准正态分布随机矩阵。 这几个函数的调用格式相似，如果这个函数的参数只是一个，那么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= 0 0 0 0 0 0 0 0 0 （2）建立一个3 2的零矩阵。 zeros(3,2) （3）设A为2 3矩阵，则可以用zeros(size(A)建立一个与矩阵A同样大小的零矩阵。 A=1 2 3;4 5 6; %产生一个23阶矩阵A zeros(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】将101125等25个数填入一个5行5列的表格中，使其每行、每列及对角线的和均为565。 一个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=vander(1:4) A = 1 1 1 1 8 4 2 1 27 9 3 1 64 16 4 1,3希尔伯特矩阵 希尔伯特（Hilbert）矩阵是一种数学变换矩阵，它的每个元素hij = 1/(i + j 1)。在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）矩阵除第一行第一列外，其他每个元素都与左上角的元素相同。生成托普利兹矩阵的函数是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是一个多项式的系数向量，高次幂系数排在前，低次幂排在后。例如，为了求多项式的x3 7x + 6的伴随矩阵，可使用如下命令： p=1,0,-7,6; compan(p) ans= 0 7 -6 1 0 0 0 1 0,6帕斯卡矩阵 我们知道，二次项(x + y)n展开后的系数随n的增大组成一个三角形表，称为杨辉三角形。由杨辉三角形表组成的矩阵称为帕斯卡（Pascal）矩阵，它的元素p1j = 1，p i1 = 1，p ij = p i 1, j 1 + p i 1, 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 + y4,4.2 矩阵分析,4.2.1 矩阵结构变换 1对角阵 只有对角线上有非0元素的矩阵称为对角矩阵，对角线上的元素都为1的对角矩阵称为单位矩阵。 （1）提取矩阵的对角线元素 设A为m n矩阵，函数diag(A)用于提取矩阵A主对角线元素，产生一个具有min(m,n)个元素的列向量。例如： A=1,2,3;4,5,6 A= 1 2 3 4 5 6 D=diag(A) D= 1 5 diag(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 = 0 0 0 0 1 0 0 0 0 2 0 0 0 0 3 0 省略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中，转置运算符是单撇号()。 A=1:3; 1:3; 1:3 B=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 0 rot90(A,4) ans = 9 37 38 -2 31 8 0 84 5,5矩阵的翻转 矩阵的翻转分左右翻转和上下翻转。对矩阵实施左右翻转是将原矩阵的第一列和最后一列调换，第二列和倒数第二列调换依此类推。 对矩阵A实施左右翻转的函数是fliplr(A)。 对矩阵A实施上下翻转的函数是flipud(A)。,4.2.2 矩阵求值 1方阵的行列式值 把一个方阵看作一个行列式，并对其按行列式的规则求值，这个值就称为矩阵所对应的行列式的值。在MATLAB中，求方阵A所对应的行列式的值的函数是det(A)。例如： A=1:3;2:-1:0;12,5,9 A= 1 2 3 2 1 0 12 5 9 B=det(A) B= -33,2矩阵的秩与迹 （1）矩阵的秩 矩阵线性无关的行数和列数称为矩阵的秩。 rank(A) A=1,2,3;2,5,6;3,2,5; r=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 范数下的条件数。 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个特征向量的线性组合表示。,例如： 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,4.3 线性方程组求解,4.3.1 矩阵求逆及线性代数方程组求解 1矩阵求逆 若方阵A、B满足等式：A B = B A = I ( I为单位阵 )，称A为B的逆矩阵 函数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.2000 A*B ans = 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)*b x= -0.3333 0.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=Ab1 y=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分解、Cholesky分解以及奇异分解等。 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分解： L,U,P=lu(A),clear A=2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4; b=13,-9,6,0; tic x2=Ab; %用左除运算求解 toc tic x1=inv(A)*b; %用求逆运算求解 toc tic L,U=lu(A); %LU分解 x3=U(Lb); %用LU分解求解 toc x1=x1 x2=x2 x3=x3,其中tic和toc两个函数配合使用用于计算程序的执行时间，tic记录当前时间，toc记录或显示使用tic函数以来所花费的时间。运行结果说明，x1、x2、x3的值相同，通过LU分解求值所化运行时间最少。,4.4.2 矩阵的QR分解 对矩阵X进行QR分解，就是把X分解为一个正交矩阵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.5774 R = -2.4495 -1.2247 -1.6330 -5.3072 0 -2.1213 2.8284 0.7071 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为一个阶数为q = p 1的上三角阵，且满足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*R ans = 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等，这些函数名都在上述内部函数名之后缀以m，并规定输入参数A必须是方阵。例如： A=4,2;3,6; B=sqrtm(A) B= 1.9171 0.4652 0.6978 2.3823 2通用矩阵函数funm funm(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.0000 expm(A) ans = 5.4366 -2.7183 2.7183 -0.0000 注意：funm函数可以用于exp函数和log函数，但求矩阵的平方根只能用sqrtm。,4.6 稀疏矩阵的处理,稀疏矩阵中具有大量的零元素，而仅含极少量的非零元素。 4.6.1 矩阵存储方式 1完全存储方式: 矩阵的全部元素按列存储。以前讲到的矩阵的存储方式都是按这个方式存储的，此存储方式对稀疏矩阵也适用。例如，不论是mn阶普通的还是稀疏的实矩阵均需要mn个存储单元，而复矩阵还要翻倍。在这种方式下，矩阵中的全部零元素也必须输入。 2稀疏存储方式: 仅存储矩阵所有的非零元素的值及其位置，即行号和列号。,4.6.2 矩阵的稀疏存储方式 1将完全存储方式转化为稀疏存储方式 函数A = sparse(S)将矩阵S转化为稀疏存储方式的矩阵A。当矩阵S是稀疏存储方式时，则函数调用相当于A = S。,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)：生成一个m n的所有元素都是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元素的下标和元素。这里产生的u、v、S可作为sparse(u,v,S)的参数。 full(A)：返回和稀疏存储矩阵A对应的完全存储方式矩阵。,2产生稀疏存储矩阵 B = spconvert(A) 其中A为一个m 3或m 4的矩阵，其每行表示一个非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) B = (3,1) -1 (2,2) 1