第7章-线性代数在后续课程中的应用举例..ppt

第七章第七章 线性代数在后续课程中的应用线性代数在后续课程中的应用举例举例¡7.1  电路电路¡7.2  信号与系统课程信号与系统课程¡7.3  数字信号处理数字信号处理¡7.4  静力学静力学¡7.5  运动学运动学¡7.6  测量学测量学¡7.7  文献管理文献管理¡7.8  经济管理经济管理7.1 电路例例7.1  电阻电路的计算电阻电路的计算 如图7-1的电路已知： 2, =4, =12, =4, =12, =4, =2，设电压源 =10v，求 ， ， 图7-1 应用实例7.1的电路图 解：用回路电流法进行建模选回路如图，设每个网孔 按图可列出各的回路电流分别为 ， 和 据根基尔霍夫定律，任何回路中诸元件上的电压之和等于零回路的电压方程为： 写成矩阵 形式为：  也可把参数代入，直接列出数字方程 简写成 A*I b* us 其中I[ ； ； ]今 =10，解此矩阵方程，即可得问题的解。

因此，可列出MATLAB程序ea701  A=[18,-12,0; -12,28,-12; 0,-12,18]; b=[1;0;0]; us=10; U=rref([A,b*us]) 程序运行的结果为： 意味着 7.2  信号与系统课程信号与系统课程例例7.3  线性时不变系统的零输入响线性时不变系统的零输入响   描述n阶线性时不变（LTI）连续系统的微分方程为 n≥m 已知y及其各阶导数的初始值为y(0)，y(1)(0)，…，y(n-1)(0)，求系统的零输入响应解： 建模 当LTI系统的输入为零时，其零输入响应为微分方程的齐次解（即令微分方程等号右端为0），其形式为（设特征根均为单根） 其中p1，p2，…，pn是特征方程a1n+a2n-1+…+ an+ an+1 =0的根，它们可用roots(a)语句求得。

各系数C1，…，Cn由y及其各阶导数的初始值来确定对此有 C1+ C2+…+Cn = y0 y0 = y(0) p1C1+ p2C2+…+ pnCn=Dy0 (Dy0表示y的导数的初始值y(1)(0)) ………………………………… 写成矩阵形式为  即 V·C = Y0 其解为 C =V \ Y0 式中 V为范德蒙矩阵，在MATLAB的特殊矩阵库中有vander函数可直接生成  MATLAB程序ea703.m a=input('输入分母系数向量a=[a1，a2，...]= '); n=length(a)-1; Y0=input(‘输入初始条件向量 Y0=[y0，Dy0， D2y0，...]= '); p=roots(a);V=rot90(vander(p));c= V\Y0'; dt=input('dt='); tf=input('tf= ') t=0:dt:tf; y=zeros(1，length(t)); for k=1:n y= y+c(k)*exp(p(k)*t);end plot(t，y)，grid  程序运行结果 用这个通用程序来解一个三阶系统，运行此程序并输入 a=[3，5，7，1]; dt=0.2; tf=8; 而Y0取 [1，0，0];[0，1，0];[0，0，1] 三种情况，用hold on语句使三次运行生成的图形画在一幅图上，得到图7-3。

图7-3 三阶系统的零输入响应7.3  数字信号处理数字信号处理   例例7.5  数字滤波器系统函数数字滤波器系统函数[12] 数字滤波器的网络结构图实际上也是一种信号流图它的特点在于所有的相加节点都限定为双输入相加器；另外，数字滤波器器件有一个迟延一个节拍的运算，它也是一个线性算子，它的标注符号为z 1根据这样的结构图，也可以用类似于例7.4的方法，求它的输入输出之间的传递函数，在数字信号处理中称为系统函数 图7-5 某数字滤波器结构图 图7-5表示了某个数字滤波器的结构图，现在要求出它的系统函数，即输出y与输入u之比先在它的三个中间节点上标注信号的名称 , , ，以便对每个节点列写方程由于迟延算子z 1不是数，要用符号代替，所以取q z 1，按照图示情况，可以写出： 写成矩阵形式为经过移项后，系统函数W可以写成：  现在可以列写计算系统函数的MATLAB程序ea705， syms q Q(1,2)q; Q(2,3)=3/8*q1/4; Q(3,1)=1; Q(3,3)=0; P=[2;1/4;0] W=inv(eye(3)Q)*P程序运行的结果为W = [16/(83*q^22*q)2*q/(83*q^22*q)] [ 2*(3*q2)/(83*q^22*q)2/(83*q^22*q)][16/(83*q^22*q)2*q/(83*q^22*q)] 以yx3作为输出的系统函数，故再键入 pretty(W(3)) 整理后得到 7.4  静力学静力学 例7.6 求双杆系统的支撑反力 图7-6 两杆系统的受力图（左）和分立体受力图（右） 两杆系统受力图见图7-6，设已知：G1=200; G2=100; L1= 2; L2=sqrt(2);theta1 =30*pi/180; theta2 = 45*pi/180; 求所示杆系的支撑反力Na,Nb,Nc  解解：◆画出杆1和杆2的分离体图，其中Na,Nb,Nc都用其x,y方向的分量Nax,Nay, Nbx,Nby, Ncx,Ncy表示，于是可列出方程如下： 对杆件1， ΣX=0 Nax + Ncx = 0 ΣY=0 Nay + Ncy - G1 = 0; ΣM=0 Ncy*L1*cos(theta1)-Ncx*L1*sin(theta1)-G1*L1/2*cos(theta1)=0; 对杆件2 ΣX=0 Nbx - Ncx = 0; ΣY=0 Nby - Ncy - G2 = 0; ΣM=0 Ncy*L2*cos(theta2)+ Ncx*L2*sin(theta2)+G2*L2/2*cos(theta2)=0; 这是一组包含六个未知数Nax, Nay, Nbx, Nby, Ncx, Ncy的六个线性代数方程,用手工解时要寻找简化的方法,用了MATLAB工具，就可直接列出矩阵方程 ，（其中X为列矩阵 ,用矩阵除法来解。

◆ MATLAB程序ea706 G1=200; G2=100; L1= 2; L2 = sqrt (2) ;% 给原始参数赋值 theta1 = 30*pi/180; theta2 = 45*pi/180;% 将度化为弧度 % 设则按此次序，系数矩阵A,B可写成下式 A=[1,0,0,0,1,0;0,1,0,0,0,1;0,0,0,0,-sin(theta1), cos(theta1);... 0,0,1,0,-1,0;0,0,0,1,0,-1;0,0,0,0, sin(theta2),cos(theta2)]; B=[0;G1;G1/2*cos(theta1);0;G2;-G2/2*cos(theta2)]; X = A\B % 用左除求解线性方程组 ◆程序运行的结果为： 即 解毕Nax Nay Nbx Nby Ncx Ncy 95.0962 154.9038 -95.0962145.0962 -95.096245.0962 7.5  运动学运动学   例例7.7  刚体空间运动学刚体空间运动学           以飞行器为例，它在空中可以围绕三个轴旋转。

假如它在向北飞行，机头正对北方，则它围绕铅垂轴的旋转角称为偏航角（Yaw），它描述了飞机左右的偏转，用u表示；围绕翼展轴的旋转角称为倾斜角（Pitch），它描述了飞机俯仰姿态，用v表示；围绕机身轴的旋转角称为滚动角（Roll），用w表示；u,v和w三个变量统称为欧拉角，它们完全描述了飞机的姿态 MATLAB中有一个演示程序quatdemo.m，专门演示这几个姿态角所造成的飞机状态键入： >>quatdemo 图7-7 飞行器姿态角演示 屏幕上将出现图7-7的画面左方为飞行器在三维空间中的模型，其中红色的是飞行器右上方为三个姿态角u,v,w的设定标尺和显示窗，右下方为在地面坐标系中的另外三个姿态角：方位角、俯仰角和倾侧角左下方还有【静态】和【动态】两个复选钮，我们只介绍【静态】，读者可自行试用【动态】进行演示 用键入参数或移动标尺的方法分别给u,v,w赋值并回车后，就可以得出相应的飞行器姿态，同时出现一根蓝色的线表示合成旋转的转轴。

 用数值来讨论这个程序的实现方法先把飞行器的三维图像用N个顶点的三维坐标描述，写成一个3N的数据矩阵G其顶点次序要适当安排，使得用plot3命令时按顶点连线能绘制出飞行器的外观例如以下的程序（ea707前半部分）即可画出一个最简单的飞行器立体图，如图7-8所示 图7-8 用数据集G画出的飞行器 Gw=[4,3,0;4,3,0;0,7,0;4,3,0]'; Gt=[0,3,0;0,3,3;0,2,0;0,3,0]'; G=[Gw,Gt] plot3(Gw(1,:),Gw(2,:),Gw(3,:),'r'),hold on plot3(Gt(1,:),Gt(2,:),Gt(3,:),'g'), axis equal 运行此程序得出整个飞行器外形的数据集为 （7.1）  飞行器围绕各个轴的旋转的结果，表现为各个顶点坐标发生变化，也就是G的变化。

只要把三种姿态的变换矩阵Y,P和R乘以图形数据矩阵G即可其中 （7.2） （7.3） （7.4） 单独变化某个姿态角所生成的图形由G1 Y*G，，G2 P*G，，G3 R*G算出，如果同时变化三个姿态角，则最后的图像数据成为Gf Y*P*R*G Q*G这里假定转动的次序为：先滚动R，再倾斜P，最后偏航Y，由于矩阵乘法不服从交换律，转动次序不同时结果也不同最后变换矩阵成为 （7.5）  此程序ea707后半部分如下： syms u w v Y=[cos(u),sin(u),0;sin(u cos(u),0;0,0,1)] R=[1,0,0;0,cos(w),sin(w);0,sin(w),cos(w)] P=[cos(v),0,sin(v);0,1,0;sin(v),0,cos(v)] Q=Y*P*R 这里采用了符号运算工具箱以得到普遍的公式，当设定了u,v,w的具体数值后，比如uπ/4,v0,wπ/3，要求出Q矩阵的数字结果时，要用subs（代换）命令如下： Asubs(Q,{u,v,w},{pi/4,0,pi/3}) 运行结果为： 我们知道， 二维变换矩阵的行列式代表的是两个向量组成的平行四边形的面积，三维变换矩阵的行列式代表的是三个向量组成的平行六面体的体积。

不难算出，det(Y)1，det(R)1，det(P)1，det(Q)1当然也有det(A)1说明这几个变换都不会改变图形对象的体积这是描述刚体运动的一个必须遵守的原则不仅不允许改变体积，而且不允许改变形状，后者要求其变换的特征值必须等于1如果特征值是复数，则它们的绝对值必须为1至于对特征向量的要求，读者可从二维情况中的A5得到启发和推论因为求特征值的MATLAB函数eig还不能用于符号函数，检验上述结论只能用数值矩阵A键入[p,lamda]eig(A) 得到p 0.7701 0.1833  0.4122i 0.1833 0.4122i 0.3190 0.6702 0.6702 0.5525 0.1315  0.5745i 0.13150.5745ilamda 1.0000 0 0 0 0.2803  0.9599i 0 0 0 0.2803  0.9599i 再键入abs(lamda)，得到 ans 1.0000 .0 . 0 0 1.0000 . 0 0 .01.0000 说明所有的特征值绝对值均为1。

为了计算各特征向量的范数，可以键入p'*p，得到ans 1.0000 0.0000  0.0000i 0.00000.0000i0.0000  0.0000i 1.0000 0.00000.0000i0.0000  0.0000i 0.0000  0.0000i 1.0000 可见其特征向量的范数也均为1 飞行器在空中的运动应该由六个自由度表征，其中三个是转动，三个是平移，要完整地描述这些运动，三维空间和33的变换矩阵是肯定不够的所以就需要研究三维以上的线性空间和线性变换问题就本题来看，研究的对象不是整个飞行器，而是飞行器外形上的N个顶点这些顶点用在三维空间中的三个坐标来描述，也就是3N数据集G3考虑了平移运动后，如同二维情况那样，也必须要用齐次坐标系，即把三维空间扩展一维，在四维空间来建立数据组和44变换矩阵即数据集G4成为4N矩阵，其最下面一行的元素都取值为1 (7.6) 而平移变换矩阵M和旋转变换矩阵Y,R,P成为 (7.7) (7.8) 其中 , , 为在三个轴 , , 方向上的平移距离。

这种方法在机器人运动学研究中也适用7.6  测量学测量学   例例7.8  用坐标测量仪测定圆形工件用坐标测量仪测定圆形工件          精密三坐标测量仪可以测量物体表面上任何一点的三维坐标，人们根据这些点的坐标就能推算出物体的其他特征尺寸比如为了测量一个圆形截面的半径，要在x-y平面内测量其圆周上n个点的坐标( , )(i=1,…n)，然后拟合出其最小二乘圆的半径 设圆周方程为： 其中为圆心的坐标，为半径对整理上述方程，得到 （7.9） 其中 ，因而 ，求出就可求出r。

 用n个测量点坐标( , )代入，得到 （7.10） 这是关于三个未知数的n个线性方程，所以是一个超定问题解出就可得知这个最小二乘圆的圆心坐标和半径r的值： 举一个数字例，设测量了圆周上7个点，其x，y坐标如下：x = -3.000 -2.000 -1.000 0 1.000 2.000 3.000y= 3.03 3.90 4.35 4.50 4.40 4.02 3.26试求此工件的直径、及圆心坐标解：7点应该有7个行向的方程，其结构相同，只是数据不同可以把数据写成列向量，然后用元素群运算一次列出所有的7个方程其程序为： x=[-3:3]'; y=[3.03,3.90,4.35,4.50,4.40,4.02,3.26]'; A=[2*x, 2*y, ones(size(x))] B=x.^2+y.^2 c=inv(A'*A)*A'*B, r=sqrt(c(3)+c(1)^2+c(2)^2)程序运行的最后结果为： c = 0.1018 工件圆心的x坐标 0.4996 工件圆心的y坐标 15.7533 r = 4.0017工件的半径r 为了使更加清楚，把运行此程序前四行得出A和B的结果写成方程： 可以看出把MATLAB的元素群运算和矩阵运算相结合，可以把复杂的计算变得何等简明。

7.7  文献管理文献管理   例例7.9  情报检索模型情报检索模型           假如数据库中包括了n个文件，而搜索所用的关键词有m个如果关键词按字母顺序排列，我们就可以把数据库表示为mn的矩阵A每个文件用矩阵的列表示A的第j列的第一个元素是一个数，它表示第一个关键词出现的相对频率，第二个元素表示第二个关键词出现的相对频率，…依此类推用于搜索的关键词清单用 空间的向量x表示如果关键词清单中第i个关键词在搜索行中出现，则x的第i个元素就赋值1，否则就赋值0为了进行搜索，只要把 乘以x 举例来说，假如我们的数据库包含有以下的书名： B1应用线性代数 B2，初等线性代数 B3初等线性代数及其应用 B4，线性代数及其应用 B5线性代数及应用 B6，矩阵代数及应用 B7矩阵理论 而搜索的6个关键词组成的集由以下的拼音字母次序排列： 初等，代数，矩阵，理论，线性，应用 因为这些关键词在书名中出现最多只有1次，所以其相对频率数不是0就是1。

当第i个关键词出现在第j本书名上时，元素A（i，j）就等于1，否则就等于0这样我们的数据库矩阵就可用下表表示：关键词 书 B1B2B3B4B5B6B7初等 0110000代数 1111110矩阵 0000011理论 0000001线性 1111100应用 1011110 假如输入的关键词是“应用，线性，代数”，则数据库矩阵和搜索向量为： 搜索结果可以表示为两者的乘积 ，于是          y的各个分量就表示各书与搜索向量匹配的程度因为    3，说明四本书 , , , 必然包含所有三个关键词这四本书就被认为具有最高的匹配度，因而在搜索的结果中把这几本书排在最前面 本例把线性变换的概念进一步扩展它不一定是在具体的几何空间内进行变量的变换（或映射），在本例中，它是从“关键词”子空间变换为“文献目录”的子空间在信号处理中，可以把时域信号变换为频域的频谱，也属于线性变换 7.8  经济管理经济管理   例例7.10  宏观经济模型宏观经济模型 假定一个国家的经济可以分解为n个部门，这些部门都有生产产品或服务的独立功能。

设单列n元向量x是这些n个部门的产出，组成在 空间的产出向量考虑它有外部的需求，要向外提供产品设d为外部需求向量，它也是 空间的单列向量，表示这n个部门以外的非经济部门对该国经济各部门的需求，比如政府消费、出口和战略储备等等 在各经济部门进行满足外部需求的生产时，它们也必须增加内部的相互需求，这种各部门的内部交叉需求非常复杂Leontiff提出的问题是，为了满足外部的最终需求向量d，各生产部门的实际产出x应该是多少，这对于经济计划的制订当然很有价值因为 x={内部需求}{外部需求d} Leontiff的输入输出模型中的一个基本假定是：对于每个部门，存在着一个在 空间单位消耗列向量 ，它表示第i个部门每产出一个单位（比如100万美金）产品，需要消耗其他部门产出的数量把这n个 ，并列起来，它可以构成一个nn的系数矩阵，可称为内部需求矩阵V由于要向外部提供产品，故内部需求矩阵各列向量元素的和必然小于1（上例中要求它等于1） 举一个最简单的例子，假如国民经济由三个部门组成，它们是制造业、农业和服务业。

它们的单位消耗列向量如下表向下列部门购买 每单位输出的输入消耗 制造业 农业服务业制造业 0.50.40.2农业 0.20.350.15服务业 0.150.10.3 如果制造业产出了100个单位的产品，有50个单位会被自己消耗，20个单位被农业消耗，而被服务业消耗的是15个单位，用算式表示为 这就是内部消耗的计算方法，把几个部门都算上，可以写出 内部需求  其中 于是总的需求方程可以写成为：               x – Vx = d               (I – V )x=d  从而可用MATLAB语句写出其解的表示式              x = inv( I –V )*d   用一个数字来试算一下，设外部需求为  用以下程序ea710解出x V=[0.5,0.4,0.2;0.2,0.35,0.15;0.15,0.1,0.3] d = [30;20;10] x = inv( eye(3) - V ) * d 程序运行的结果为 这个结果是合理的，因为实际产出应该比外部需求大得多，以应付内部的消耗。

我们常常说，某某外部需求可以拉动国民经济增长多少个百分点，就是从这样的模型中得出的 内部需求矩阵V要满足一些基本要求，一般各列的列向元素总和必须小于1，否则这个部门就将入不敷出而亏损，但我们仍可能求出上述方程的解当所有的列向量都出现列向元素总和大于1的情况时，解x中会出现负值，因而是庸解请注意分析请注意分析其实际意义其实际意义。