《线性代数Matlab教案第7章课件.ppt》由会员分享，可在线阅读，更多相关《线性代数Matlab教案第7章课件.ppt（52页珍藏版）》请在三一办公上搜索。
1、7.1 电路,例7.1 电阻电路的计算 如图7-1的电路。已知：2,=4,=12,=4,=12,=4,=2，设电压源=10v，求，。图7-1 应用实例7.1的电路图,解：用回路电流法进行建模。选回路如图，设每个网孔。按图可列出各的回路电流分别为，和。据根基尔霍夫定律，任何回路中诸元件上的电压之和等于零回路的电压方程为：写成矩阵 形式为：,也可把参数代入，直接列出数字方程 简写成 A*Ib*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*
2、us)程序运行的结果为：意味着,7.2 信号与系统课程,例7.3 线性时不变系统的零输入响 描述n阶线性时不变（LTI）连续系统的微分方程为 nm 已知y及其各阶导数的初始值为y(0)，y(1)(0)，y(n-1)(0)，求系统的零输入响应。,解：建模 当LTI系统的输入为零时，其零输入响应为微分方程的齐次解（即令微分方程等号右端为0），其形式为（设特征根均为单根）其中p1，p2，pn是特征方程a1n+a2n-1+an+an+1=0的根，它们可用roots(a)语句求得。各系数C1，Cn由y及其各阶导数的初始值来确定。对此有,C1+C2+Cn=y0 y0=y(0)p1C1+p2C2+pnCn=
3、Dy0(Dy0表示y的导数的初始值y(1)(0)写成矩阵形式为,即 VC=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=VY0;dt=input(dt=);tf=input(tf=)t=0:dt:tf;y=zeros(1，length(t);for k=1:n y=y+c(k)*ex
4、p(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
5、的方法，求它的输入输出之间的传递函数，在数字信号处理中称为系统函数。,图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*q1/4;Q(3,1)=1;Q(3,3)=0;P=2;1/4;0 W=inv(eye(3)Q)*P程序运行的结果为
6、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)整理后得到,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都
7、用其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的六个线性代数方程,用手工解时要寻找简
8、化的方法,用了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
9、);0;G2;-G2/2*cos(theta2);X=AB%用左除求解线性方程组,程序运行的结果为：即 解毕。,7.5 运动学,例7.7 刚体空间运动学 以飞行器为例，它在空中可以围绕三个轴旋转。假如它在向北飞行，机头正对北方，则它围绕铅垂轴的旋转角称为偏航角（Yaw），它描述了飞机左右的偏转，用u表示；围绕翼展轴的旋转角称为倾斜角（Pitch），它描述了飞机俯仰姿态，用v表示；围绕机身轴的旋转角称为滚动角（Roll），用w表示；u,v和w三个变量统称为欧拉角，它们完全描述了飞机的姿态。,MATLAB中有一个演示程序quatdemo.m，专门演示这几个姿态角所造成的飞机状态。键入：quatde
10、mo 图7-7 飞行器姿态角演示,屏幕上将出现图7-7的画面。左方为飞行器在三维空间中的模型，其中红色的是飞行器。右上方为三个姿态角u,v,w的设定标尺和显示窗，右下方为在地面坐标系中的另外三个姿态角：方位角、俯仰角和倾侧角。左下方还有【静态】和【动态】两个复选钮，我们只介绍【静态】，读者可自行试用【动态】进行演示。用键入参数或移动标尺的方法分别给u,v,w赋值并回车后，就可以得出相应的飞行器姿态，同时出现一根蓝色的线表示合成旋转的转轴。,用数值来讨论这个程序的实现方法。先把飞行器的三维图像用N个顶点的三维坐标描述，写成一个3N的数据矩阵G。其顶点次序要适当安排，使得用plot3命令时按顶点连
11、线能绘制出飞行器的外观。例如以下的程序（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乘以
12、图形数据矩阵G即可。其中（7.2）（7.3）（7.4）,单独变化某个姿态角所生成的图形由G1Y*G，G2P*G，G3R*G算出，如果同时变化三个姿态角，则最后的图像数据成为GfY*P*R*GQ*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)
13、,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)运行结果为：,我们知道，二维变换矩阵的行列式代表的是两个向量组成的平行四边形的面积，三维变换矩阵的行列式代表的是三个向量组成的平行六面体的体积。不难算出，det(Y)1，det(R)1，det(P)1，det(Q)1。当然也有det(A)1。说明这几个变换都不会改变图形对象的体积。这是描述刚体运动的一个必须遵守的原则。不仅不允许改变体积，而且不允许改变形状
14、，后者要求其变换的特征值必须等于1。如果特征值是复数，则它们的绝对值必须为1。至于对特征向量的要求，读者可从二维情况中的A5得到启发和推论。因为求特征值的MATLAB函数eig还不能用于符号函数，检验上述结论只能用数值矩阵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)，得到 an
15、s 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.0000i0.0000 0.0000i 1.0000 可见其特征向量的范数也均为1。,飞行器在空中的运动应该由六个自由度表征，其中三个是转动，三个是平移，要完整地描述这些运动，三维空间和33的变换矩阵是肯定不够的。所以就需要研究三维以上的线性空间和线性变换问题。就本题来看，研究的对象不是整
16、个飞行器，而是飞行器外形上的N个顶点。这些顶点用在三维空间中的三个坐标来描述，也就是3N数据集G3。考虑了平移运动后，如同二维情况那样，也必须要用齐次坐标系，即把三维空间扩展一维，在四维空间来建立数据组和44变换矩阵。即数据集G4成为4N矩阵，其最下面一行的元素都取值为1。,(7.6)而平移变换矩阵M和旋转变换矩阵Y,R,P成为(7.7)(7.8)其中,为在三个轴,方向上的平移距离。这种方法在机器人运动学研究中也适用。,7.6 测量学,例7.8 用坐标测量仪测定圆形工件 精密三坐标测量仪可以测量物体表面上任何一点的三维坐标，人们根据这些点的坐标就能推算出物体的其他特征尺寸。比如为了测量一个圆形
17、截面的半径，要在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点应
18、该有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的元素群运算和矩阵运算相结合，
19、可以把复杂的计算变得何等简明。,7.7 文献管理,例7.9 情报检索模型 假如数据库中包括了n个文件，而搜索所用的关键词有m个。如果关键词按字母顺序排列，我们就可以把数据库表示为mn的矩阵A。每个文件用矩阵的列表示。A的第j列的第一个元素是一个数，它表示第一个关键词出现的相对频率，第二个元素表示第二个关键词出现的相对频率，依此类推。用于搜索的关键词清单用 空间的向量x表示。如果关键词清单中第i个关键词在搜索行中出现，则x的第i个元素就赋值1，否则就赋值0。为了进行搜索，只要把 乘以x。,举例来说，假如我们的数据库包含有以下的书名：B1。应用线性代数 B2，初等线性代数 B3。初等线性代数及其应
20、用 B4，线性代数及其应用 B5。线性代数及应用 B6，矩阵代数及应用 B7。矩阵理论 而搜索的6个关键词组成的集由以下的拼音字母次序排列：初等，代数，矩阵，理论，线性，应用,因为这些关键词在书名中出现最多只有1次，所以其相对频率数不是0就是1。当第i个关键词出现在第j本书名上时，元素A（i，j）就等于1，否则就等于0。这样我们的数据库矩阵就可用下表表示：,假如输入的关键词是“应用，线性，代数”，则数据库矩阵和搜索向量为：搜索结果可以表示为两者的乘积，于是,y的各个分量就表示各书与搜索向量匹配的程度。因为 3，说明四本书,必然包含所有三个关键词。这四本书就被认为具有最高的匹配度，因而在搜索的结
21、果中把这几本书排在最前面。本例把线性变换的概念进一步扩展。它不一定是在具体的几何空间内进行变量的变换（或映射），在本例中，它是从“关键词”子空间变换为“文献目录”的子空间。在信号处理中，可以把时域信号变换为频域的频谱，也属于线性变换。,7.8 经济管理,例7.10 宏观经济模型 假定一个国家的经济可以分解为n个部门，这些部门都有生产产品或服务的独立功能。设单列n元向量x是这些n个部门的产出，组成在 空间的产出向量。考虑它有外部的需求，要向外提供产品。设d为外部需求向量，它也是 空间的单列向量，表示这n个部门以外的非经济部门对该国经济各部门的需求，比如政府消费、出口和战略储备等等。,在各经济部门
22、进行满足外部需求的生产时，它们也必须增加内部的相互需求，这种各部门的内部交叉需求非常复杂。Leontiff提出的问题是，为了满足外部的最终需求向量d，各生产部门的实际产出x应该是多少，这对于经济计划的制订当然很有价值。因为 x=内部需求外部需求d Leontiff的输入输出模型中的一个基本假定是：对于每个部门，存在着一个在 空间单位消耗列向量，它表示第i个部门每产出一个单位（比如100万美金）产品，需要消耗其他部门产出的数量。把这n个，并列起来，它可以构成一个nn的系数矩阵，可称为内部需求矩阵V。由于要向外部提供产品，故内部需求矩阵各列向量元素的和必然小于1（上例中要求它等于1）。,举一个最简
23、单的例子，假如国民经济由三个部门组成，它们是制造业、农业和服务业。它们的单位消耗列向量如下表。,如果制造业产出了100个单位的产品，有50个单位会被自己消耗，20个单位被农业消耗，而被服务业消耗的是15个单位，用算式表示为 这就是内部消耗的计算方法，把几个部门都算上，可以写出 内部需求,其中 于是总的需求方程可以写成为：xVx=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中会出现负值，因而是庸解。请注意分析其实际意义。,
