线性代数实践教师班第三讲.ppt

线性代数实践（第三讲）,第7章 矩阵运算法解方程,7.1 矩阵运算的规则,在MATLAB入门中已讲过的，不再重复。 由于其乘法不符合交换律，有些公式不能乱用； 单列向量与单行向量的左右两种乘法要加区别，而且往往有特别的用途。 例如向量长度（范数）的计算； 例如二维坐标网格的生成； X=ones(21,1)*-10:10,Y= -10:10*ones(1,21) 矩阵的乘幂An, eA和(I-A) -1的级数展开，都要求A是方阵。,矩阵乘法不满足交换律,有许多我们习惯的公式，其中隐含地包含了交换律，这些公式在矩阵运算中也不能直接使用。比如： 正确的做法是展开时不交换次序,平面上网格坐标系的产生,用列矩阵乘行矩阵生成网格坐标,这两个矩阵都是21行21列的，都有441个元素，如何快捷地输入呢？这时可以用到列乘行的乘法运算。可用下面的语句： h10:10; lhlength(h) % 输入均分行向量 % 用全么列乘均分行生成X Xones(lh,1)*h % 用均分列乘全么行生成Y Yh*ones(1,lh),7.2 初等变换乘子矩阵的生成,行交换E1gen(n,i,j)：使n行矩阵中的第i,j两行交换 function E=E1gen(n,i,j) n=size(A); E=eye(n); E(i,i)=0; E(j,j)=0; E(i,j)=1; E(j,i)=1; 乘子矩阵E2gen (n,i,k),使n行矩阵中的第i行乘以k function E=E2gen(n,i,k) n=size(A); E=eye(n); E(i,i)=k; E3gen(n,i,j,c)使n行矩阵中的第i行乘以k加到第j行上 function E=E3gen(n,i,j,k) n=size(A); E=eye(n); E(j,i)=k;,初等变换乘子矩阵示例,E=E1gen(8,4,6) E2=E2gen(8,4,6) E3=E3gen(8,4,6,5) 例如E3=E3gen(3,1,3,4),例7.2.4 求消元所需的乘子矩阵,要消去下列矩阵的A(2,1),求乘子矩阵E3 在第二行加以第一行乘A(2,1)/A(1,1)3， 故令B E3gen（A,1,2,3）,行阶梯生成等价于矩阵左乘,因此，整个行阶梯形式U的生成过程，可以看作把原矩阵左乘以一系列的初等变换矩阵E1和E3。把这些初等矩阵的连乘积写成Ex，设其逆为L: 从而有 L*UA （7.10） 就是说，A可以分解为一个准下三角矩阵L和一个上三角（即行阶梯）矩阵U的乘积。MATLAB提供了三角分解的函数lu，它的调用方法是： L,Ulu(A),lu分解是求行阶梯的一个方法,用lu函数求出的U实际上就是A的行阶梯形式（不是简化行阶梯形式）。所以，求简化行阶梯形式用rref函数，而求行阶梯形式可以用lu 函数。不过，它和我们用消元运算所得U的数据不一定相同，尽管得出的阶次和阶梯形状相同。但因为行阶梯形式可以有无数种，用不同步骤算出的结果也不同。只有变成简化行阶梯形式，才能进行比较，看它是不是惟一的。,7.3 行列式的定义和计算,两种定义方法： 1。按全排列求和定义， 其中tj为第j种排列的逆序数。,行列式第2种定义方法,2。按解的分母项，从低阶到高阶用归纳法定义 二阶： 三阶：,两种定义方法的比较,第一种定义的两个数学难点全排列和逆序数，是绝大多数工科学生一生不会用的。 第二种定义方法自然地得出了行列式按行（或按列）展开的公式。美国教材都用第二种定义方法，成电教材（全国精品课程）也用这种方法。 两种方法都不能用来计算，因为其计算效率都极低，2525矩阵要算上万年。 第8章将指出，行列式的几何意义是面积或体积，可否从这方面探索，因为它的用途很单一，就是判断奇异性，连正负号都不必关心。,行列式的计算方法,计算行列式的最好方法还是行阶梯法，可以利用lu分解 L,Ulu(A) 把A分解为一个准下三角矩阵L和一个上三角矩阵U的乘积。因为det(L)1，所以U和A的行列式相等。 det(A) det(U) 而三角矩阵U的det(U)很好求。只要把U的主对角线元素连乘就可得到它的行列式。 此法所需的乘法次数仅为定义1法的10 -23,行列式计算实例7.3.1,程序如下 l,ulu(A), du diag(u) Dprod(du) 结果为,du 10 4.8 10.625 9.4824 1.2349 D 5.9720e003 5972,7.4 矩阵的秩和矩阵求逆,按定义，矩阵的秩是矩阵A中行列式不等于零的最高阶子式的阶次。是用以衡量联立方程中有效方程数目的指数。 按照定义来计算矩阵的秩，可能遇到的问题也是子矩阵的数量很大，每个矩阵的行列式计算又非常麻烦，其计算量也将是不可接受的天文数字。 计算矩阵的秩的最好方法仍然是行阶梯法，如第6章所述，行阶梯化简后非全为零的行数，就是该矩阵的秩。用MATLAB函数rrank(A)可以检验A的秩，rank函数对A是否是方阵没有要求，即可以有mn。,矩阵求逆,对于nn方阵A，当rn时，称A是满秩的，若rn，必有det(A)0，称A是欠秩的或奇异的。奇异矩阵不可以求逆。 矩阵求逆的最简单方法也是行阶梯化简，其方法是设定一个由A和I组成的增广矩阵CA,I，求C的简化行阶梯形式UCrref(A,I)，得出UCI,V。V就显示出这个逆矩阵的内容。,例7.6 求逆矩阵示例,求A的逆阵 解：程序ag706。 A3,0,3,6;5,1,1,5;3,1,4,9;1,3,4,4; CA,eye(4) U0Crref(C) VU0C(:,5:8),程序运行结果,右边四列就是其逆阵： 矩阵求逆命令：V=inv(A),用inv函数求逆,求A的逆阵 程序ag707为： A=-16,-4,-6;15,-3,9;18,0,9, V=inv(A) 运行结果： Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.042030e-018.,条件数衡量奇异程度的量,在用数值方法计算矩阵的逆时，由于计算中的误差，人们不大可能得到理想的零合理想的全零行，所以矩阵是否奇异，并不是那么绝对的。 为了评价矩阵接近奇异的程度，采用了条件数（Condition Number）作为常用的衡量指标。它永远大于1。其数值愈接近于1，计算误差愈小； MATLAB中，条件数用cond(A)计算，它达到104以上时，求逆的误差就可能相当可观。像现在，条件数达到1016（注：条件数是逆条件数RCOND的倒数），结果是根本不能用的。,7.5 用矩阵除法解线性方程,如果mn，则线性代数方程 Axb （7.21） 中的A是方阵，设det(A)0，则它的逆阵存在。将上式左右同乘以inv(A) ，由于inv(A) *AI，得到 xinv(A)*b （7.23） MATLAB创立了矩阵除法的概念，因为 inv(A)相当于将A放到分母上去，所以可以把上式写成 xA b （7.24） 就称为左除，因为inv(A)是乘在b的左方。,左除解线性方程的扩展,左除的功能远远超过了矩阵求逆函数inv，inv(A)函数要求A必须是方阵，所以（7.23）式只能用来解适定方程，而（7.24）式并不要求A为方阵，在A是mn阶且mn（欠定）时，它只要求A与b的行数相等且A的秩为m。所以（7.24）式也可以用来解欠定方程，在下例中可以看出。此外，运算符还能用来解超定方程，,左除解欠定方程,例7.8 用矩阵算法解例6.5.1 A3,4,3,2,1;0,6,0,3,3;4,3,4,2,2; 1,1,1,0,1;2,6,2,1,3; b 2; 3; 2; 0; 1; x=Ab 得到x=inf，无解。改用行阶梯方法找有效行。 左除要求的是系数矩阵的行数与秩相同， BA,b,r=rank(B),UB,iprref(B); U0UB(1:r,1:5); dUB(1:3,6); xU0d,本例运行结果,r=3, 及 它是此欠定方程的一个特解。,7.6.1 网络的矩阵分割和连接,在电路设计中，经常要把复杂的电路分割为局部电路，每一个电路都用一个网络黑盒子来表示。黑盒子的输入为u1，i1，输出为u2，i2，其输入输出关系用矩阵A来表示（如图7.1所示）： A是22矩阵，称为该局部电路的传输矩阵,两个网络的串联,两个串接的子网络。第一个子网络包含电阻R1，第二个子网络包含电阻R2，列出第一个子网络的电路方程为： 由 得矩阵方程,两个网络的串联（续）,由第二网络： 写成矩阵方程为： 整个电路的传输矩阵为两者的乘积,7.6.2 用逆阵进行保密编译码,在英文中有一种对消息进行保密的措施，就是把英文字母用一个整数来表示。然后传送这组整数。这种方法是很容易根据数字出现的频率来破译，例如出现频率特别高的数字，很可能对应于字母E。 可以用乘以矩阵A的方法来进一步加密。假如A是一个行列式等于1的整数矩阵，则A1的元素也必定是整数。而经过这样变换过的消息，同样两个字母对应的数字不同，所以就较难破译。 接收方只要将这个消息乘以A1就可以复原。,7.6.3 减肥配方的实现,设脱脂牛奶的用量为x1个单位（100g），大豆面粉的用量为x2个单位，乳清的用量为x3个单位，表中的三个营养成分列向量为： 使这个合成的营养与剑桥配方的要求相等，得到,7.6.4 弹性梁的柔度矩阵,设简支梁如图7.3所示，在梁的三个位置分别施加力f1，f2和f3后，在该处产生的综合变形为图示的y1，y2和y3，通常称为挠度。根据虎克定律，在材料未失去弹性的范围内，力与它引起的变形呈线性关系，可以写出： 矩阵中的元素d为单位力f引起的挠度，它愈大，表明这个梁愈柔软。,数字实例,设柔度矩阵 （1）在1，2，3处施加的力为30，50和20试求出其挠度。 （2）要在3处产生0.4挠度，其他两处为零，求应加的力。 程序ag764 D0.001*5,2,1;2,4,3;1,3,6 % 输入柔度矩阵 f30;50;20, yD*f （排齐）% 给定力，求挠度 y10;0;0.4 % 给定挠度， Kinv(D), f1K*y1 % 求刚度矩阵,求力,梁的刚度矩阵计算,柔度矩阵的逆就是刚度矩阵K，KD 1， 其中,7.6.5 网络和图,图为1，2，3，4四个城市之间的空运航线，用有向图表示。则该图可以用下列航路矩阵表示： 经过一次转机（也就是坐两次航班）能到达的城市，可以由邻接矩阵的平方A2A12来求得。,第8章 用向量空间解方程组,8.1 向量和向量空间 二维空间R2中的向量用两个沿列向的元素表示 u=2;4; v=3;-1; plot(2,3,4,1,x);hold on % 若用中的子程序drawvec， drawvec(u);hold on drawvec(v,g);hold off,二维向量张成的空间,平面上的任何一点w1;w2是不是一定能用u和v的线性组合来实现？即是不是一定能找到一组常数c1,c2，使得 c1,c2取所有可能的值，得到的w的集合就是u和v张成的子空间，在所给的u和v下，它是一个平面。 若u和v两个向量的各元素成简单的比例关系，合成的向量只能在一根直线上，不可能张成整个二维平面。这种情况下，称这两个向量u和v是线性相关的。,2三维空间中的向量,若v1,v2和v3都是三维空间的列向量。可以用空间坐标中的三个点，或从坐标原点引向这三点的箭头来表示。用矩阵代数表示如下 如果三个基本向量之间线性无关，那么它们的线性组合可以覆盖（张成）整个三维空间。如果三个向量共面，即相关，就不能张成三维空间。判断三个向量的线性相关性，可用行列式。,三维空间向量的相关性,即看三向量并列所得矩阵的行列式 det(A)=0 相关 det(A)0 不相关 行列式的几何意义:在二维是两个向量组成的平行四边形面积，在三维是三个向量组成的平行六面体的体积。,行列式的几何意义,二维 三维 det(A)=右图平行六面体的体积,n维向量的相关性,在进入三维以上的空间时，已经没有可与面积、体积直接相当的概念可用了，所以采用了秩的概念。如果A的行列式为零，也就是它的秩r小于n时，说明这n个向量是线性相关的。 秩的概念也概括了面积存在（r2）和体积存在（r3）的意义，因此，它是更高度的抽象。,8.2 向量空间和基向量,若r个向量是线性无关的，则它们的线性组合的全体V就构成了r维空间Rr 。如果它不是空集，则V称为向量空间。生成V的r个线性无关的向量v称为基向量或基（Basis）。 当rn时，给定的n个向量就是一组基。如果rn，那就要在n个向量中选出r个线性无关的向量。用秩的概念还无法判定哪些向量是线性无关的，这时又要藉助于把矩阵简化为阶梯形式的方法。,例8.2 求四个五维向量的子空间,这四个向量组成的矩 阵如右，对它进行行 阶梯简化。程序为： A4,5,4,1;0,3,0,1;2,1,2,0;5,4,5,3;1,4,1,1 U0,iprref(A) 得到 ip=1,2,4 其三个枢轴列对应的就是 三个线性无关的列向量。,三个向量空间位置演示程序,三维空间中，为了观察三个向量的空间关系，ATLAST手册还提供了一个演示程序viewsubspaces(u,v,w)，它用蓝色直线显示向量u，同时用红色显示v和w所组张成的平行四边形平面，画在同一张立体图上。例如： u=-1;1;8;v=5;-4;7;w=-3;1;-5; viewsubspaces(u,v,w),grid on 三个向量的起点都是xyz0的原点。要看清其几何意义，还是需要一定的空间想象力。,三个向量的空间关系,例8.3 w是否在v1,v2,v3的空间内,设 w是否能由v1,v2,v3的线性组合构成的问题，取决于线性方程组 解的存在性。 v1=7;-4;-2;9; v2=-4;5;-1;-7; v3=9;4;4;-7; w=-9;7;1;-4; v=v1,v2,v3; c=vw % 把基向量组成矩阵v求解 也可以按det(v)是否为零进行判别,8.3 向量的内积和正交性,在三维空间中，x和y两个向量的内积定义为x,yx1y1x2y2x3y3。m维情况可以写成 这是一个标量。向量x与自己求内积： 得到的是其各分量的平方和，其平方根就等于向量的长度（或模、或范数norm）。,内积的几何意义,在平面情况，两向量的内积除以它们的长度是它们夹角的余弦，可以利用下图证明。 根据余弦定律， 最后得到 此结果可推广到高维空间，只是被抽象化了：,例8.4 基向量长度规一化和夹角,例8.4 求例8.3中的单位基向量v01,v02,v03，并分别求它们之间的夹角。 解：解题的程序为ag822： v10=v1/norm(v1), v20=v2/norm(v2), v30=v3/norm(v1), theta12=acos(v1*v2)/(norm(v1)*norm(v2) theta13=acos(v1*v3)/(norm(v1)*norm(v3) theta23=acos(v3*v2)/(norm(v3)*norm(v2) ),正交基向量的生成,两向量x，y正交的条件是它们的内积为零。 给出向量求正交基常用施密特算法，ATLAST手册中给出了相应的程序gschmidt。调用时键入Q,R=gschmidt(v)，Q就是单位正交基向量e。 MATLAB中不用施密特算法，而用更好的算法编成了正交分解子程序qr.m，它将v分解为Q和R两个矩阵的乘积。调用方法为： Q,Rqr(v) Q就是mm单位正交矩阵。,基向量正交化的schmidt公式,得到qi(i1,2,k)后，再把它们除以norm(qi)，就可归一化为单位向量ek。,基向量正交化的schmidt子程序,function Q,R=gschmidt(V) m,n=size(V); R=zeros(n); R(1,1)=norm(V(:,1); Q(:,1)=V(:,1)/R(1,1); for k=2:n R(1:k1,k)=Q(:,1:k1)*V(:,k); Q(:,k)=V(:,k)Q(:,1:k1)*R(1:k1,k); R(k,k)=norm(Q(:,k); Q(:,k)=Q(:,k)/R(k,k); end,求单位正交基向量的例,例8.5 对于例8.3的数据，求其规范化正交基向量e1,e2,en。 解：程序为 V7,4,9;4,5,4;2,1,4;9,7,7 Q,Rqr(v) % 或 Q,Rgschmidt(v) eQ(:,1:3) 得到：,8.4 齐次方程Ax=0的解空间,设有m个方程和n个变量，A的秩是r，则经过行简化后得到的行阶梯矩阵U的有r个枢轴元素，非枢轴元素有nr个。因此该方程的全解将等于Axb的一个特解加上其齐次方程Ax0的通解。本节将从向量空间的视点来讨论它的解，因为通解是nr阶的无穷的集合，所以要研究解所张成的向量空间。 Ax0意味着这些解x的集合经过矩阵A变换后都映射到像空间的零点，所以英文把此解所张成的空间称为Null Space，直译为零空间。我国的通用译名为解空间或基础解系，我们觉得用齐次解空间较为准确。,齐次方程Ax=0解空间的例,例8.6 试求下列系数矩阵的齐次解空间： 解：输入A，并求出它的简化行阶梯形式，键入u0,iprref(A)，得到 ip1,3,齐次解空间的例（续）,其通解可以看成三个向量的线性组合 这个式子就表示了一个三维的向量空间，在这个空间中所有的向量都能使Ax0。所以它被称为齐次解空间或零空间。,求齐次解空间的子程序,这样齐次解空间的m is系数矩阵N可以用下面的程序来自动完成： functin N=nulspace(A) m,n=size(A); U0,ip=rref(A) is=1:n; is(ip)=; N(ip,:)=-U0(1:rank(A),is); N(is,:)=eye(n-rank(A) MATLAB中的子程序为N=null(A).,计算例题8.7,系数矩阵A如右， 求Ax0的通解。 解：程序ag842 先输入A，再键入 vnulbasis(A) %或 v=null(A,r) r表示用有理分式的 基向量 得到都是三个分量并列， v=v1,v2,v3,8.5 解超定方程的思路,有时用向量空间的方法可以更为简捷地推导公式，超定方程的解就是一个例子。 既然我们已讨论了适定方程组和不定方程组的求解方法，自然会提出如何解超定方程组的问题。 工程问题都可以允许方程有误差，把一组解代入方程后，每个方程都有误差；要找误差在一定意义下的总和为最小的解。在这样的思路引导下，就产生了超定方程求解的方法。,误差线性方程组的建立,引入误差向量e。 eAxb 写出其完全的矩阵形式如下 问题是，找到解x，使e的长度或范数为最小。,从向量空间的视点分析,研究例6.1的超定方程组（d）： 改写成 简写为 选择不同的x1和x2将得到不同的合成向量A*xx1v1x2v2q ，q必定处于v1和v2张成的平面之内。而方程中的b则一般不会在这个平面内，,本例的向量空间图,这时最近似的解就应该是该平面上与b点最近的点所对应的坐标A*xhat。它应该是b点向v1和v2张成的平面的投影。所以和b的连线应该和v1和v2张成的平面垂直，也就是说必须分别与v1和v2正交。如图8.6所示。,最小二乘解的公式推导,A*xhat和b的连线向量应该是这两个向量之差，即 ，它与v1和v2正交的要求可以分别表示为： 和 综合在一起可以写成： 最后得到公式,最小二乘解的数字例,例8. 求题6.1(d)方程组的最小二乘解。 解：MATLAB程序ag808如下： A=1,1;1,1; 1,2, b=1;3;3 xhat=inv(A*A)*A*b e=A*xhatb, norm(e) 运行此程序，得到,MATLAB中超定方程的解,在MATLAB中，把运算(ATA) -1AT单独编成一个子程序，称为pinv函数。求最小二乘解的公式可以写成 xpinv(A)*b，与适定方程的解xinv(A)*b非常相似，只是pinv函数并不要求A是方阵。 最小二乘解也可用运算符表示，这就把欠定方程、适定方程和超定方程用统一的运算格式： xAb MATLAB会自动根据系数矩阵A的行数m和列数n，来判断采用哪个方法和程序。不过对于欠定方程，这个式子只给出了一个特解，没给通解。,数字实例8.9：实验数据处理,例8.9 设在某一实验中，给某元件加1,2,3,4,5v电压，测得的电流为0.2339，0.3812，0.5759，0.8153，0.9742ma。求此元件的电阻。 解：设直线的方程为y c(1)x c(2)，待定的系数是c(1)，c(2)。将上述数据分别代入x，y，把这五个方程联立，用矩阵表述：,实验数据处理实例,写成 datax *c(1) ones(N,1) * c(2) datay 其中datax，datay都是5行数据列向量，这是5个一次代数方程，含两个未知数，方程数超过未知数的数目，是一个超定方程，解的程序如下： A datax , ones(N,1); B datay; c A B,8.6.1 价格平衡模型,单位消耗列向量vi表示第i个部门每产出一个单位产品中，本部门和其他各个部门消耗的百分比。 于是总的价格平衡方程可以写成为： ( I V ) p =0