1、7.1 电路例例7.1 电阻电路的计算电阻电路的计算 如图7-1的电路。已知：2,=4,=12,=4,=12,=4,=2，设电压源 =10v，求 ，。图7-1 应用实例7.1的电路图1R2R3R4R5R6R7RSu3i4u7u 解：用回路电流法进行建模。选回路如图，设每个网孔。按图可列出各的回路电流分别为 ，和 。据根基尔霍夫定律，任何回路中诸元件上的电压之和等于零回路的电压方程为：写成矩阵 形式为：123333455556701000abscRRRRiRRRRRiuRRRRi aibici1233334555567()()0()0absabcbcRRR iR iuR iRRR iR iR i
2、RRR i 也可把参数代入，直接列出数字方程 简写成 A*I b*us 其中I ；。今 =10，解此矩阵方程，即可得问题的解。因此，可列出MATLAB程序ea701 1812011228120.012180absciiui aibicisu A=18,-12,0;-12,28,-12;0,-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
3、线性时不变系统的零输入响线性时不变系统的零输入响 描述n阶线性时不变（LTI）连续系统的微分方程为 nm 已知y及其各阶导数的初始值为y(0)，y(1)(0)，y(n-1)(0)，求系统的零输入响应。,dddddddddd111121ubtubtubyatyatyatyammmmnnnnn解：建模 当LTI系统的输入为零时，其零输入响应为微分方程的齐次解（即令微分方程等号右端为0），其形式为（设特征根均为单根）其中p1，p2，pn是特征方程a1n+a2n-1+an+an+1=0的根，它们可用roots(a)语句求得。各系数C1，Cn由y及其各阶导数的初始值来确定。对此有1212()eeenp
4、tp tp tny tCCC C1+C2+Cn=y0 y0=y(0)p1C1+p2C2+pnCn=Dy0 (Dy0表示y的导数的初始值y(1)(0)写成矩阵形式为 111111220DnnnnnnpCpCpCy0100211121121DD111yyyCCCppppppnnnnnnn 即 VC=Y0 其解为 C=V Y0 式中 V为范德蒙矩阵，在MATLAB的特殊矩阵库中有vander函数可直接生成。112000,;,D,DnnC CCyyyTT0CY1121121111nnnnnppppppV MATLAB程序ea703.m a=input(输入分母系数向量a=a1，a2，.=);n=len
5、gth(a)-1;Y0=input(输入初始条件向量 Y0=y0，Dy0，D2y0，.=);p=roots(a);V=rot90(vander(p);c=VY0;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语句使三次运行生成的图形画在一幅图上，得到图
6、7-3。图7-3 三阶系统的零输入响应7.3 数字信号处理数字信号处理 例例7.5 数字滤波器系统函数数字滤波器系统函数12 数字滤波器的网络结构图实际上也是一种信号流图。它的特点在于所有的相加节点都限定为双输入相加器；另外，数字滤波器器件有一个迟延一个节拍的运算，它也是一个线性算子，它的标注符号为z 1。根据这样的结构图，也可以用类似于例7.4的方法，求它的输入输出之间的传递函数，在数字信号处理中称为系统函数。图7-5 某数字滤波器结构图 图7-5表示了某个数字滤波器的结构图，现在要求出它的系统函数，即输出y与输入u之比。先在它的三个中间节点上标注信号的名称 ,，以便对每个节点列写方程。由于
7、迟延算子z 1不是数，要用符号代替，所以取q z 1，按照图示情况，可以写出：1x2x3x写成矩阵形式为经过移项后，系统函数W可以写成：1223312311844xqxuxqxuxx112233002311008440100qxxxqxxx 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*q
8、22*q)2/(83*q22*q)16/(83*q22*q)2*q/(83*q22*q)以yx3作为输出的系统函数，故再键入 pretty(W(3)整理后得到 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都用
9、其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的六个线性代数方程,用手工解时要寻找简化的
10、方法,用了MATLAB工具，就可直接列出矩阵方程 ，（其中X为列矩阵 ,用矩阵除法来解。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(thet
11、a2);B=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 Nb
12、x 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，专门演示这几个姿态角所造
13、成的飞机状态。键入：quatdemo 图7-7 飞行器姿态角演示 屏幕上将出现图7-7的画面。左方为飞行器在三维空间中的模型，其中红色的是飞行器。右上方为三个姿态角u,v,w的设定标尺和显示窗，右下方为在地面坐标系中的另外三个姿态角：方位角、俯仰角和倾侧角。左下方还有【静态】和【动态】两个复选钮，我们只介绍【静态】，读者可自行试用【动态】进行演示。用键入参数或移动标尺的方法分别给u,v,w赋值并回车后，就可以得出相应的飞行器姿态，同时出现一根蓝色的线表示合成旋转的转轴。用数值来讨论这个程序的实现方法。先把飞行器的三维图像用N个顶点的三维坐标描述，写成一个3N的数据矩阵G。其顶点次序要适当安排，
14、使得用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 2 30 0 0 0 0 3
15、 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*R*G Q*G。这里假定转动的次序为：先滚动R，再倾斜P，最后偏航Y，
16、由于矩阵乘法不服从交换律，转动次序不同时结果也不同。最后变换矩阵成为 （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()*cos(w)vv 此程序ea707后半部分如下：syms u w v Y=cos(u),sin(u),0;sin(u cos(u),0;0,0,1
17、)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.7071 0.3536 0.6124 0 0.8660 0.5000 A 我们知道，二维变换矩阵的行列式代表的是两个向量组成的平行四边形的面积，三维变换矩阵的行列
18、式代表的是三个向量组成的平行六面体的体积。不难算出，det(Y)1，det(R)1，det(P)1，det(Q)1。当然也有det(A)1。说明这几个变换都不会改变图形对象的体积。这是描述刚体运动的一个必须遵守的原则。不仅不允许改变体积，而且不允许改变形状，后者要求其变换的特征值必须等于1。如果特征值是复数，则它们的绝对值必须为1。至于对特征向量的要求，读者可从二维情况中的A5得到启发和推论。因为求特征值的MATLAB函数eig还不能用于符号函数，检验上述结论只能用数值矩阵A。键入p,lamdaeig(A)得到p 0.7701 0.1833 0.4122i 0.1833 0.4122i 0.3
19、190 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 0.0000 0.0000i 0.00000.0000i0.0000 0.0000i 1.0000 0.00000.0000i0.0000 0.0000i 0.0000 0.0
20、000i 1.0000 可见其特征向量的范数也均为1。飞行器在空中的运动应该由六个自由度表征，其中三个是转动，三个是平移，要完整地描述这些运动，三维空间和33的变换矩阵是肯定不够的。所以就需要研究三维以上的线性空间和线性变换问题。就本题来看，研究的对象不是整个飞行器，而是飞行器外形上的N个顶点。这些顶点用在三维空间中的三个坐标来描述，也就是3N数据集G3。考虑了平移运动后，如同二维情况那样，也必须要用齐次坐标系，即把三维空间扩展一维，在四维空间来建立数据组和44变换矩阵。即数据集G4成为4N矩阵，其最下面一行的元素都取值为1。(7.6)而平移变换矩阵M和旋转变换矩阵Y,R,P成为 (7.7)(
21、7.8)其中 ,为在三个轴 ,方向上的平移距离。这种方法在机器人运动学研究中也适用。34GG=ones(1,N)123100010,0010001ccc4M0000001333444Y,(R,P)Y,(R,P)1x2x3x1c2c3c7.6 测量学测量学 例例7.8 用坐标测量仪测定圆形工件用坐标测量仪测定圆形工件 精密三坐标测量仪可以测量物体表面上任何一点的三维坐标，人们根据这些点的坐标就能推算出物体的其他特征尺寸。比如为了测量一个圆形截面的半径，要在x-y平面内测量其圆周上n个点的坐标(,)(i=1,n)，然后拟合出其最小二乘圆的半径。ixiy 设圆周方程为：其中为圆心的坐标，为半径。对整
22、理上述方程，得到 （7.9）其中 ，因而 ，求出就可求出r。22212()()xcycr22222121212322()22xcycr ccxcyc cxy 222312crcc22312rccc 用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.262211111222222232222122
23、1221nnnnxyxycxyxyccxyxyAc=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(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 为了使更加清楚，
24、把运行此程序前四行得出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.0000 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 情报检索模型情
25、报检索模型 假如数据库中包括了n个文件，而搜索所用的关键词有m个。如果关键词按字母顺序排列，我们就可以把数据库表示为mn的矩阵A。每个文件用矩阵的列表示。A的第j列的第一个元素是一个数，它表示第一个关键词出现的相对频率，第二个元素表示第二个关键词出现的相对频率，依此类推。用于搜索的关键词清单用 空间的向量x表示。如果关键词清单中第i个关键词在搜索行中出现，则x的第i个元素就赋值1，否则就赋值0。为了进行搜索，只要把 乘以x。TAmR 举例来说，假如我们的数据库包含有以下的书名：B1。应用线性代数 B2，初等线性代数 B3。初等线性代数及其应用 B4，线性代数及其应用 B5。线性代数及应用 B6
26、，矩阵代数及应用 B7。矩阵理论 而搜索的6个关键词组成的集由以下的拼音字母次序排列：初等，代数，矩阵，理论，线性，应用 因为这些关键词在书名中出现最多只有1次，所以其相对频率数不是0就是1。当第i个关键词出现在第j本书名上时，元素A（i，j）就等于1，否则就等于0。这样我们的数据库矩阵就可用下表表示：关键词 书 B1B2B3B4B5B6B7初等 0110000代数 1111110矩阵 0000011理论 0000001线性 1111100应用 1011110 假如输入的关键词是“应用，线性，代数”，则数据库矩阵和搜索向量为：搜索结果可以表示为两者的乘积 ，于是 0 1 1 0 0 0 0 1
27、 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=A x 30 21 30 30 31 1 0 0 1 21 0 0 1 1 0 0 0 TyA x y的各个分量就表示各书与搜索向量匹配的程度。因为 3，说明四本书 ,必然包含所有三个关键词。这四本书就被认为具有最高的匹配度，因而在搜索的结果中把这几本书排在最前面。本例把线性变换的概念进一步扩展。它不一定是在
28、具体的几何空间内进行变量的变换（或映射），在本例中，它是从“关键词”子空间变换为“文献目录”的子空间。在信号处理中，可以把时域信号变换为频域的频谱，也属于线性变换。1y3y4y5y1B3B4B5B7.8 经济管理经济管理 例例7.10 宏观经济模型宏观经济模型 假定一个国家的经济可以分解为n个部门，这些部门都有生产产品或服务的独立功能。设单列n元向量x是这些n个部门的产出，组成在 空间的产出向量。考虑它有外部的需求，要向外提供产品。设d为外部需求向量，它也是 空间的单列向量，表示这n个部门以外的非经济部门对该国经济各部门的需求，比如政府消费、出口和战略储备等等。nRnR 在各经济部门进行满足外
29、部需求的生产时，它们也必须增加内部的相互需求，这种各部门的内部交叉需求非常复杂。Leontiff提出的问题是，为了满足外部的最终需求向量d，各生产部门的实际产出x应该是多少，这对于经济计划的制订当然很有价值。因为 x=内部需求外部需求d Leontiff的输入输出模型中的一个基本假定是：对于每个部门，存在着一个在 空间单位消耗列向量 ，它表示第i个部门每产出一个单位（比如100万美金）产品，需要消耗其他部门产出的数量。把这n个 ，并列起来，它可以构成一个nn的系数矩阵，可称为内部需求矩阵V。由于要向外部提供产品，故内部需求矩阵各列向量元素的和必然小于1（上例中要求它等于1）。nRiviv 举一
30、个最简单的例子，假如国民经济由三个部门组成，它们是制造业、农业和服务业。它们的单位消耗列向量如下表。向下列部门购买 每单位输出的输入消耗 制造业 农业服务业制造业 0.50.40.2农业 0.20.350.15服务业 0.150.10.3 如果制造业产出了100个单位的产品，有50个单位会被自己消耗，20个单位被农业消耗，而被服务业消耗的是15个单位，用算式表示为 这就是内部消耗的计算方法，把几个部门都算上，可以写出 内部需求 0.5501001000.2200.15151v11 12 13 312323,xxvxvxvv v vxx V x 其中 于是总的需求方程可以写成为：x Vx=d (
31、I V)x=d 从而可用MATLAB语句写出其解的表示式 x=inv(I V)*d 用一个数字来试算一下，设外部需求为 0.50.40.20.20.35 0.150.150.10.3V302010d 用以下程序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 程序运行的结果为 这个结果是合理的，因为实际产出应该比外部需求大得多，以应付内部的消耗。我们常常说，某某外部需求可以拉动国民经济增长多少个百分点，就是从这样的模型中得出的。160.4563 94.4867 62.1673x 内部需求矩阵V要满足一些基本要求，一般各列的列向元素总和必须小于1，否则这个部门就将入不敷出而亏损，但我们仍可能求出上述方程的解。当所有的列向量都出现列向元素总和大于1的情况时，解x中会出现负值，因而是庸解。请注意分析请注意分析其实际意义。其实际意义。
