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

1、第七章第七章 线性代数在后续课程中的应用线性代数在后续课程中的应用举例举例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的电路图1R2R3R4R5R6R7RSu3i4u7u 解：用回路电流法进行建模。选回路如图，设每个网孔。 按图可列出各的回路电流分别为 ，

2、 和 。据根基尔霍夫定律，任何回路中诸元件上的电压之和等于零回路的电压方程为： 写成矩阵 形式为： 123333455556701000abscRRRRiRRRRRiuRRRRi aibici1233334555567()()0()0absabcbcRRR iR iuR iRRR iR iR iRRR i 也可把参数代入，直接列出数字方程 简写成 A*I b* us 其中I ； ； 。今 =10，解此矩阵方程，即可得问题的解。 因此，可列出MATLAB程序ea701 1812011228120.012180absciiui aibicisu A=18,-12,0; -12,28,-12; 0,

3、-12,18; b=1;0;0; us=10; U=rref(A,b*us) 程序运行的结果为： 意味着 1.0000 0 0 0.9259 0 1.0000 0 0.5556 0 0 1.0000 0.3704U0 0.9259 0.5556 0.3704abciiiI7.2 信号与系统课程信号与系统课程例例7.3 线性时不变系统的零输入响线性时不变系统的零输入响 描述n阶线性时不变（LTI）连续系统的微分方程为 nm 已知y及其各阶导数的初始值为y(0)，y(1)(0)，y(n-1)(0)，求系统的零输入响应。,dddddddddd111121ubtubtubyatyatyatyammmm

4、nnnnn解： 建模 当LTI系统的输入为零时，其零输入响应为微分方程的齐次解（即令微分方程等号右端为0），其形式为（设特征根均为单根） 其中p1，p2，pn是特征方程a1n+a2n-1+ an+ an+1 =0的根，它们可用roots(a)语句求得。各系数C1，Cn由y及其各阶导数的初始值来确定。对此有1212( )eeenp tp tp tny tCCC C1+ C2+Cn = y0 y0 = y(0) p1C1+ p2C2+ pnCn=Dy0 (Dy0表示y的导数的初始值y(1)(0) 写成矩阵形式为 111111220DnnnnnnpCpCpCy0100211121121DD111yy

5、yCCCppppppnnnnnnn 即 VC = Y0 其解为 C =V Y0 式中 V为范德蒙矩阵，在MATLAB的特殊矩阵库中有vander函数可直接生成。 112000, ;,D ,DnnC CCyyyTT0CY1121121111nnnnnppppppV 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= VY0; dt=input(dt=); tf=input(tf= )

6、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 数字滤波器的网络结构图实际上也是一种信号流图。它的特点在于所有的相加节点都限

7、定为双输入相加器；另外，数字滤波器器件有一个迟延一个节拍的运算，它也是一个线性算子，它的标注符号为z 1。根据这样的结构图，也可以用类似于例7.4的方法，求它的输入输出之间的传递函数，在数字信号处理中称为系统函数。 图7-5 某数字滤波器结构图 图7-5表示了某个数字滤波器的结构图，现在要求出它的系统函数，即输出y与输入u之比。先在它的三个中间节点上标注信号的名称 , , ，以便对每个节点列写方程。由于迟延算子z 1不是数，要用符号代替，所以取q z 1，按照图示情况，可以写出： 1x2x3x写成矩阵形式为经过移项后，系统函数W可以写成： 1223312311844xqxuxqxuxx1122

8、33002311008440100qxxxqxxx xux = Qx-PuW = x/u = inv(I -Q)*P 现在可以列写计算系统函数的MATLAB程序ea705， syms q Q(1,2)q; Q(2,3)=3/8*q1/4; Q(3,1)=1; Q(3,3)=0; P=2;1/4;0 W=inv(eye(3)Q)*P程序运行的结果为W = 16/(83*q22*q)2*q/(83*q22*q) 2*(3*q2)/(83*q22*q)2/(83*q22*q)16/(83*q22*q)2*q/(83*q22*q) 以yx3作为输出的系统函数，故再键入 pretty(W(3) 整理后得

9、到 1222116 288(3)8 321.541.54yqqzWuqqqqzz 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

10、= 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为列矩

11、阵 ,用矩阵除法来解。 MATLAB程序ea706 G1=200; G2=100; L1= 2; L2 = sqrt (2) ;% 给原始参数赋值 theta1 = 30*pi/180; theta2 = 45*pi/180;% 将度化为弧度AX = BNax;Nay; Nbx; Nby; Ncx; Ncy % 设则按此次序，系数矩阵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

12、=0;G1;G1/2*cos(theta1);0;G2;-G2/2*cos(theta2); X = AB % 用左除求解线性方程组 程序运行的结果为： 即 解毕。 1 0 0 0 1.0000 0 0 1 0 0 0 1.0000 0 0 0 0 -0.5000 0.8660 0 0 1 0 -1.0000 0 0 0 A 0 95.0962 200.0000 154.903 86.6025, 0 0 1 0 -1.0000 100.0000 0 0 0 0 0.7071 0.7071 -35.3553BX8 -95.0962 145.0962 -95.0962 45.0962Nax Nay

13、 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，专门演示这几

14、个姿态角所造成的飞机状态。键入： quatdemo 图7-7 飞行器姿态角演示 屏幕上将出现图7-7的画面。左方为飞行器在三维空间中的模型，其中红色的是飞行器。右上方为三个姿态角u,v,w的设定标尺和显示窗，右下方为在地面坐标系中的另外三个姿态角：方位角、俯仰角和倾侧角。左下方还有【静态】和【动态】两个复选钮，我们只介绍【静态】，读者可自行试用【动态】进行演示。 用键入参数或移动标尺的方法分别给u,v,w赋值并回车后，就可以得出相应的飞行器姿态，同时出现一根蓝色的线表示合成旋转的转轴。 用数值来讨论这个程序的实现方法。先把飞行器的三维图像用N个顶点的三维坐标描述，写成一个3N的数据矩阵G。其顶

15、点次序要适当安排，使得用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） 4 4 0 4 0 0 0 03 3 7 3 3 3

16、2 30 0 0 0 0 3 0 0 G 飞行器围绕各个轴的旋转的结果，表现为各个顶点坐标发生变化，也就是G的变化。只要把三种姿态的变换矩阵Y,P和R乘以图形数据矩阵G即可。其中 （7.2） （7.3） （7.4） cos( ) sin( ) 0sin( ) cos( ) 0 0 0 1uuuu Y 1 0 0 0 cos(w) sin(w) 0 sin(w) cos(w)R cos( ) 0sin( ) 0 1 0 sin( ) 0 cos( )vvvvP 单独变化某个姿态角所生成的图形由G1 Y*G，G2 P*G，G3 R*G算出，如果同时变化三个姿态角，则最后的图像数据成为Gf Y*P*

17、R*G Q*G。这里假定转动的次序为：先滚动R，再倾斜P，最后偏航Y，由于矩阵乘法不服从交换律，转动次序不同时结果也不同。最后变换矩阵成为 （7.5） cos( )*cos( ) sin( )*cos(w ) cos( )*sin( )*sin(w ) sin( )*sin(w ) cos( )*sin( )*cos(w )sin( )*cos( ) cos( )*cos(w )+sin( )*sin( )*sin(w ) cos( )*sin(w )+sin( )*sin( )*cos(w ) sin( ) uvuuvuuvQuvuuvuuvv cos( )*sin(w ), cos( )*

18、cos(w )vv 此程序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,v0,w/3，要求出Q矩阵的数字结果时，要用subs（代换）命令如下： Asubs(Q,u,v,w,pi/4,0,pi/3) 运行结果为： 0.7071 0.3536 0.61240.

19、7071 0.3536 0.6124 0 0.8660 0.5000 A 我们知道， 二维变换矩阵的行列式代表的是两个向量组成的平行四边形的面积，三维变换矩阵的行列式代表的是三个向量组成的平行六面体的体积。不难算出，det(Y)1，det(R)1，det(P)1，det(Q)1。当然也有det(A)1。说明这几个变换都不会改变图形对象的体积。这是描述刚体运动的一个必须遵守的原则。不仅不允许改变体积，而且不允许改变形状，后者要求其变换的特征值必须等于1。如果特征值是复数，则它们的绝对值必须为1。至于对特征向量的要求，读者可从二维情况中的A5得到启发和推论。因为求特征值的MATLAB函数eig还不

20、能用于符号函数，检验上述结论只能用数值矩阵A。键入p,lamdaeig(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.13150.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

21、0.0000 0.0000i 0.00000.0000i0.0000 0.0000i 1.0000 0.00000.0000i0.0000 0.0000i 0.0000 0.0000i 1.0000 可见其特征向量的范数也均为1。 飞行器在空中的运动应该由六个自由度表征，其中三个是转动，三个是平移，要完整地描述这些运动，三维空间和33的变换矩阵是肯定不够的。所以就需要研究三维以上的线性空间和线性变换问题。就本题来看，研究的对象不是整个飞行器，而是飞行器外形上的N个顶点。这些顶点用在三维空间中的三个坐标来描述，也就是3N数据集G3。考虑了平移运动后，如同二维情况那样，也必须要用齐次坐标系，即把三

22、维空间扩展一维，在四维空间来建立数据组和44变换矩阵。即数据集G4成为4N矩阵，其最下面一行的元素都取值为1。 (7.6) 而平移变换矩阵M和旋转变换矩阵Y,R,P成为 (7.7) (7.8) 其中 , , 为在三个轴 , , 方向上的平移距离。这种方法在机器人运动学研究中也适用。34GG =ones(1,N)123100010,0010001ccc4M0000001333444Y ,(R ,P )Y ,(R ,P )1x2x3x1c2c3c7.6 测量学测量学 例例7.8 用坐标测量仪测定圆形工件用坐标测量仪测定圆形工件 精密三坐标测量仪可以测量物体表面上任何一点的三维坐标，人们根据这些点的

23、坐标就能推算出物体的其他特征尺寸。比如为了测量一个圆形截面的半径，要在x-y平面内测量其圆周上n个点的坐标( , )(i=1,n)，然后拟合出其最小二乘圆的半径。ixiy 设圆周方程为： 其中为圆心的坐标，为半径。对整理上述方程，得到 （7.9） 其中 ，因而 ，求出就可求出r。 22212()()xcycr22222121212322() 22xcycr ccxcyc c xy 222312crcc22312rccc 用n个测量点坐标( , )代入，得到 （7.10） 这是关于三个未知数的n个线性方程，所以是一个超定问题。解出就可得知这个最小二乘圆的圆心坐标和半径r的值： 举一个数字例，设测

24、量了圆周上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.2622111112222222322221221221nnnnxyxycxyxyccxyxyAc = Biyix试求此工件的直径、及圆心坐标解：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

25、(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的元素群运算和矩阵运算相结合，可以把复杂的计算变得何等简明。 -6.0000 6.0600 1.0000 -4.0000 7.8000 1.0000 -2.0000 8.7000 1.0000 0 9.0000 1.0000 2.0000 8.8000 1.00

26、00 4.0000 8.0400 1.0000 6.0000 6.5200 1Ac=B123 18.1809 19.2100 19.9225 20.2500 20.3600 20.1604.0000 19.6276ccc 7.7 文献管理文献管理 例例7.9 情报检索模型情报检索模型 假如数据库中包括了n个文件，而搜索所用的关键词有m个。如果关键词按字母顺序排列，我们就可以把数据库表示为mn的矩阵A。每个文件用矩阵的列表示。A的第j列的第一个元素是一个数，它表示第一个关键词出现的相对频率，第二个元素表示第二个关键词出现的相对频率，依此类推。用于搜索的关键词清单用 空间的向量x表示。如果关键词清

27、单中第i个关键词在搜索行中出现，则x的第i个元素就赋值1，否则就赋值0。为了进行搜索，只要把 乘以x。TAmR 举例来说，假如我们的数据库包含有以下的书名： B1。应用线性代数 B2，初等线性代数 B3。初等线性代数及其应用 B4，线性代数及其应用 B5。线性代数及应用 B6，矩阵代数及应用 B7。矩阵理论 而搜索的6个关键词组成的集由以下的拼音字母次序排列： 初等，代数，矩阵，理论，线性，应用 因为这些关键词在书名中出现最多只有1次，所以其相对频率数不是0就是1。当第i个关键词出现在第j本书名上时，元素A（i，j）就等于1，否则就等于0。这样我们的数据库矩阵就可用下表表示：关键词 书 B1B

28、2B3B4B5B6B7初等 0110000代数 1111110矩阵 0000011理论 0000001线性 1111100应用 1011110 假如输入的关键词是“应用，线性，代数”，则数据库矩阵和搜索向量为： 搜索结果可以表示为两者的乘积 ，于是 0 1 1 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 1 1 1 1 1A010,0 0 01 1 0 1 1 1 1 01x 0 1 0 0 1 1 1 1 0 0 1 0 1 1 0 0 1 1 0 1 0 0 1 1 0 1 0 0 1 1 0 1 Ty = Ax 30 21 30 3

29、0 31 1 0 0 1 21 0 0 1 1 0 0 0 TyA x y的各个分量就表示各书与搜索向量匹配的程度。因为 3，说明四本书 , , , 必然包含所有三个关键词。这四本书就被认为具有最高的匹配度，因而在搜索的结果中把这几本书排在最前面。 本例把线性变换的概念进一步扩展。它不一定是在具体的几何空间内进行变量的变换（或映射），在本例中，它是从“关键词”子空间变换为“文献目录”的子空间。在信号处理中，可以把时域信号变换为频域的频谱，也属于线性变换。 1y3y4y5y1B3B4B5B7.8 经济管理经济管理 例例7.10 宏观经济模型宏观经济模型 假定一个国家的经济可以分解为n个部门，这些部门都有生产产品或服务的独立功能。设单列n元向量x是这些n个部门的产出，组成在 空间的产出向量。考虑它有外部的需求，要向外提供产品。设d为外部需求向量，它也是 空间的单列向量，表示这n个部门以外的非经济部门对该国经济各部门的需求，比如政府消费、出口和战略储备等等。nRnR 在各经济部门进行满足外部需求的生产时，它们也必须增加内部的相互需求，这种各部门的内部交叉需求非常复杂。Leontiff提出的问题是