线性代数实践(教师班第10讲).ppt

第十章 后续课矩阵建模举例,10.1 多项式插值问题,例：试求三次插值多项式，使曲线通过以下4个点:(0,3),(1,0),(2,-1),(3,6)。解：这4个点的坐标应满足三次多项式函数，代入后有：程序：A=1,0,0,0;1,1,1,1;1,2,4,8;1,3,9,27;B=3;0;-1;6;a=AB得到a=3,-2,-2,1T,即,例10.1的图形,绘图程序ezplot(3-2*t-2*t2+t3)hold on,grid onplot(0:3,3,0,-1,6,x)line(1.5,1.5,0,6)axis(-1,4,-2,8)求t=1.5处的插值函数值t1=1.5;p1=3-2*t1-2*t12+t13plot(t1,p1,o),图10.1 例10.1的插值曲线,高阶的多项式插值,在一般情况下，当给出函数在n+1个点上的值时，就可以用n次多项式对它进行插值。如果给出的点数（即方程数）大于n+1，方程组成为超定的，因而没有一个能满足方程组的解，得出的曲线将是以最小二乘意义下的误差靠近各点，于是插值就变为拟合。插值也不一定是自变量的多项式，比如圆锥曲线方程 虽然它有6个系数，若用a除以此方程两端，得到的将是有5个待定系数的方程。如果给出x-y平面上的5个点，就可以列出5个线性方程来确定这5个系数。,10.2 坐标测量仪测定的拟合,比如为了测量一个圆锥形截面的半径，可在x-y平面内测量其圆周上n个点的坐标(xi,yi)(i=1,n)，然后拟合出此截面的方程。对于每一组数据(xi,yi)，代入圆锥曲线方程，移项可得：n个点就有n个方程。其结构相同，只是数据不同。可以把数据写成列向量，然后用元素群运算一次列出所有的n个方程。,例10.2 圆锥截面方程的拟合,设测量了圆周上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 试求出此圆锥截面的方程，并求其最大最小直径。列出程序如下：x=-3:3;%把x,y赋值为列向量y=3.03,3.90,4.35,4.50,4.40,4.02,3.26;A=x.*y,y.*y,x,y,ones(size(x);B=-x.2;%列出系数矩阵A和Bc=inv(A*A)*A*B,%解超定方程得出c,程序运行结果,将程序运行生成的参数写出，应为：即截面的方程为：,例10.2 曲线,ezplot(x2+0.005*x*y+0.9214*y2-0.2228*x-0.4050*y-16.8537)画圆锥曲线hold on,plot(x,y,x)画测量点grid on,axis equal%x,y两方向取同样比例尺拟合的图形和给定的测量点，如图所示。从工程角度看，这样的测量点布局对拟合精度不大有利。,估计圆直径的方法,设圆周方程为：c1,c2为圆心的坐标，r为半径。整理上述方程，得到用n个测量点坐标(xi,yi)代入，得到这是关于三个未知数的n个线性方程，所以是一个超定问题。解出就可得知这个最小二乘圆的圆心坐标和半径r的值：,例10.2a 曲线程序及运行结果,程序ag1002ax=-3:3;%把x,y赋值为列向量y=3.03,3.90,4.35,4.50,4.40,4.02,3.26;A=2*x,2*y,ones(size(x)%求出系数矩阵A,BB=x.2+y.2c=inv(A*A)*A*B,%求超定方程的解，得出cr=sqrt(c(3)+c(1)2+c(2)2)%由c求出r程序运行的最后结果为：c=0.1018 工件圆心的x坐标0.4996 工件圆心的y坐标15.7533r=4.0017工件的半径r,按圆形估计的最小二乘图形,以上就是本题要解的超定线性方程组快速绘制此圆形的语句为：ezplot(x-0.1018)2+(y-0.4996)2-4.00172=0),10.3 天体轨道测量拟合,在适当的极坐标中，天体的位置应满足下列方程：其中为常数而是轨道的偏心率，对于椭圆，对于抛物线，而对于双曲线。对一个新发现的天体的观测得到了表中的数据，试确定其轨道的性质，并预测此天体在弧度时的位置。解：对任何一组测量数据i,ri，可列出以下对和两个变量的线性方程：,0,r,求和的程序,五个测量点共有五个方程，轨道只有两个未知数，这是一个超定问题。先用元素群运算一次列出5个方程的系数矩阵，再用超定方程解的公式。theta=0.88,1.1,1.42,1.77,2.14;测定的theta 和r 数组r=3,2.3,1.65,1.25,1.01A=ones(5,1),r.*cos(theta),B=r,方程组系数矩阵 X=inv(A*A)*A*B%解超定方程 rg=X(1)/(1-X(2)*cos(0.46)%求theta=0.46处的r运行结果为：即轨道方程为，rg=5.3098,程序ag1002a的数据和图形,%画出极坐标中轨迹的语句，ezpolar(1.4509/(1-0.8111*cos(theta),hold on,polar(theta,r,x),axis(-4,8,-5,5)%画出测定点得到图形如右。,10.4 静力学问题,静力学研究物体受力后的平衡方程。一个物体在平面上平衡，需要三个平衡方程。如果是N个物体相互作用下的平衡，那么方程的总数就是3N个。空间物体的平衡需要6个平衡方程。这个方程组通常都是线性的，所以求解就可以归结为矩阵方程。它可以避免解单个方程和单个数，只要把系数矩阵输入程序，就轻而易举地同时得出所有的解。例10.4 两杆系统受力图见图，设已知：G1=200;G2=100;L1=2;L2=sqrt(2);theta1=30*pi/180;theta2=45*pi/180;求所示杆系的支撑反力Na,Nb,Nc。,例10.4的图和方程,对杆件1，X=0 Nax+Ncx=0；Y=0 Nay+Ncy-G1=0;M=Ncy*L1*cos(theta1)-Ncx*L1*sin(theta1)-G1*L1/2*cos(theta1)=0;对杆件2X=0 Nbx-Ncx=0;Y=0 Nby-Ncy-G2=0;M=Ncy*L2*cos(theta2)+Ncx*L2*sin(theta2)+G2*L2/2*cos(theta2)=0;,例10.4的矩阵模型,这是一组包含六个未知数Nax,Nay,Nbx,Nby,Ncx,Ncy的六个线性代数方程，写成矩阵方程AX+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=AB%求解线性方程组,程序ag1004的运行结果,这样求解的方法不仅适用于全部静力学题目，而且可用于材料力学和结构力学中的静力学和静不定问题。,10.5 材料力学的静不定问题,建模：A点因各杆变形而引起的x方向位移x，y方向位移y,由几何关系,得变形方程:其中Ki为杆i的刚度系数。再加上两个力平衡方程：共有n+2个方程，其中包含n个未知力和两个待求位移x和y，方程组适定可解.。,例10.5 三杆静不定桁架,图示三杆桁架,挂一重物P=3000，设L=2m,各杆的截面积分别为A1=200e-6,A2=300e-6,A3=400e-6,材料的弹性模量E=200e9，求各杆受力的大小。解：此时应有两个平衡方程和三个位移协调方程:-N1*cos(a1)-N2-N3*cos(a3)=0;N1*sin(a1)-N3*sin(a3)=P;N1/K1=dx*cos(a1)+dy*sin(a1)N2/K2=dxN3/K3=dx*cos(a3)dy*sin(a3),程序ag1005编写,令X为列向量N1;N2;N3;dx;dy，把上述五个线性方程组列成D*X=B的矩阵形式，从而就可用MATLAB的左除语句XDB来求解。因此程序ag1005如下：%设定参数，计算杆长，刚度系数等，分行给D赋值 D(1,:)=cos(a1),cos(a2),cos(a3),0,0;D(2,:)=sin(a1),sin(a2),sin(a3),0,0;D(3,:)=1/K1,0,0,-cos(a1),-sin(a1);D(4,:)=0,1/K2,0,-cos(a2),-sin(a2);D(5,:)=0,0,1/K3,-cos(a3),-sin(a3);B=P;0;0;0;0;%给B赋值format long,X=DB%解方程组，用长格式显示结果,程序ag1005运行结果,执行此程序,显示的结果为：X=1763.40607065591(N1)591.14251029634(N2)-2995.72429657297(N3)0.00016949097(dx)0.00001970475(dy)若用普通格式显示,dy=0.0000,实际上它不是零，这可从N2不等于零推想。在MATLAB中一个矩阵中包含数值相差很大的元素时,就得采用format long来显示，或者单独显示各个元素X(1),X(5)。读者还可改变几根杆的刚度系数，看它们如何影响各杆力的分布。,10.6 二自由度机械振动,右图表示了由两个质量、弹簧及阻尼器构成的二自由度振动系统，今要在给定两个质量的初始位置和初始速度的情况下求系统的运动。,解：设x1和x2分别表示两个质量关于它们平衡位置的偏差值，则此二自由度振动系统的一般方程为：,建模的矩阵形式,可写成矩阵形式：(2)其中(3)这是一个四阶常微分方程组。给出它的初始条件（初始位置X0和初始速度Xd0）:(4)可以求出它的解，但用解析方法求解非常麻烦。如果C=0，即无阻尼情况，则系统可解耦为两种独立的振动模态，通常的书上只给出解耦简化后的解。,线性微分方程的矩阵解,基本思路是把原始方程化成四个一阶方程矩阵方程组：(5)这个方程在初始条件Y=Y0下的解为。用MATLAB表示为Y=Y0*expm(A*t)。其中expm表示把(A*t)看成矩阵来求其指数。只要把Y,Y0和A找到就行。先把方程（2）写成两个一阶矩阵方程：(6)于是，(7)因 故Y和Y0都是41单列变量；,例10.6 二维振动系统的数值例,设m1=1;m2=9;k1=4;k2=2;c1和c2可选。求在初始条件x0=1;0和 xd0=0;-1下，系统的输出x1,x2曲线。根据上面的模型可以写出程序ag1006如下。m1=1;m2=9;k1=4;k2=2;%输入各原始参数c1=input(c1=);c2=input(c2=);%输入阻尼系数x0=1;0;xd0=0;-1;tf=50;dt=0.1;%给出初始条件及时间向量M=m1,0;0,m2;K=k1+k2,-k2;-k2,k2;%构成二阶参数矩阵C=c1+c2,-c2;-c2,c2;A=zeros(2,2),eye(2);-MK,-MC;%构成四阶参数矩阵y0=x0;xd0;%四元变量的初始条件for i=1:round(tf/dt)+1%设定计算点,作循环计算 t(i)=dt*(i-1);y(:,i)=expm(A*t(i)*y0;%循环计算矩阵指数end,程序运行结果的曲线,运行此程序，输入c1=0.2,c2=0.5所得的结果见左图。从中可清楚地看到振动的两种模态。特别是x1的运动反映了两种模态的迭合。给出不同的初始条件，各模态的幅度也会变化。输入c1=0,c2=0所得的结果见右图，这时两种振动模态分别出现在两个波形中，即已经简单地解耦了,振动模态分析,键入：p,lamda=eig(A)得到 特征根lamda由两组共轭复根组成，它的虚部是振动的角频率，它的实部是它的衰减系数。所以矩阵特征根反映了两种振动模态的特征，使我们更进一步理解了“特征”两字的物理意义。,10.6 交流稳态电路的复数矩阵,稳态交流电路方程是线性的，但其系数是复数，变量是平面上的向量（也是复数），所以得到的是复数线性方程组。虽然线性代数理论只证到实数，但MATLAB中有关线性代数的函数都对实数和复数同样有效，所以都可以用在复数条件下。,例10.7 如图10-10所示电路，设Z1=-j250，Z2=250，Is=20A，负载ZL=500j 500，求负载电压。,例10.7的矩阵模型,MATLAB程序ag1007,Z1=-j*250;Z2=250;ki=0.5;Is=2+j*0;zL=500+j*500;%设定元件参数a11=1/Z1+1/Z2;a12=-1/Z2;a13=0;%设定系数矩阵Aa21=-1/Z2;a22=1/Z2-1/zL;a23=-ki;a31=1/Z1;a32=0;a33=-1;A=a11,a12,a13;a21,a22,a23;a31,a32,a33;B=1;0;0;%设定系数矩阵BX=AB*Is;%求方程解X=Ua;Ub;I1Ub=X(2)%负载电压absUb=abs(Ub)%负载电压的幅度angleUb=angle(Ub)*180/pi%负载电压的相角（度）compass(Ub);%画出Ub的向量图%compass(Is;X)%画出四个向量的图,程序运行结果为：,负载电压Ub=-2.5000e+002-7.5000e+002i幅度absUb=790.5694，相角angleUb=-108.4349Ub向量见图。键入compass(Is;X)，将 画出四个向量的图，因为Ua,Ub的数值比Is,I1大了上百倍，Is,I1在图中无法显示,10.8 线性系统零输入响应的计算,线性时不变连续系统的特性可用常微分方程表示为：求其零输入响应。解：在零输入条件u=0时，等式右端为零。系统响应的通解为：其中，p是特征方程的n个根组成的向量p1,p2,pn，其每个分量的系数Cn则由y及其各阶导数的初始条件y0,Dy0,D(n-1)y0来确定。,代入初始条件得到的矩阵方程,初始条件数应该和常数数相等，由此构成一个确定C1,Cn的线性代数方程组，写成:矩阵V由特征根向量p确定，这种矩阵称为范德蒙特矩阵。,求零输入响应程序ag1008,在MATLAB中，有生成它的函数vander(p)。它产生的矩阵与本题的矩阵排列转了90度，应该用V=rot90(vander(p)，于是即可按此思路编成MATLAB程序ag1008:a=input(输入分母系数向量a=a1,a2,.=);n=length(a)-1;Y0=input(输入初始条件 Y0=y0,Dy0,D2y0,.=);p=roots(a);V=rot90(vander(p)；%生成系数矩阵Vc=VY0;%用左除求出系数向量C%以下是计算并画出时间响应dt=input(dt=);tf=input(tf=)%输入时间数组的参数t=0:dt:tf;y=zeros(1,length(t);%生成t向量和y初始向量for k=1:n y=y+c(k)*exp(p(k)*t);end%将各分量叠加plot(t,y),grid%绘图,例10.8 数字实例,用这一个普遍程序来解一个三阶系统，设其微分方程为：初始条件为：，求零输入响应。,解：运行ag1008，按提示输入a为1,2,9,3;Y0为1,0,0;dt为0.1;tf为5；即可得到t=05秒的零输入响应曲线。再分别取Y0为0,1,0;0,0,1，用hold on语句使三次运行生成的图形叠画在一幅图上，得到右图。,10.9 计算频谱用的DFT矩阵,有限长序列x(n)(0nN-1)有N个样本值。它的傅里叶变换在频率区间(02)的N个等间隔分布的点k=2k/N(0kN-1)上也有N个样本值。这两组有限长的序列之间可以用简单的关系联系起来：其中，是一个相角为2/N的单位向量，也称为旋转因子。X(k)就称为x(n)的离散傅里叶变换(DFT)，也就是x(n)的频谱。,频谱计算的矩阵模型,式(1)可以写成矩阵的形式如下：所以信号频谱的计算，可以简单地用一个矩阵乘法来完成。信号是N1维数组，矩阵W称为DFT矩阵，它是NN维的复数阵。把矩阵乘法看作一变换，我们就可以把频谱计算看作信号从时域到频域的线性变换。N的典型值是1024。,实现频谱计算的程序(1),这个矩阵运算可用下列几条MATLAB语句来实现。n=0:1:N-1;k=n;%设定n和k的行向量WN=exp(-j*2*pi/N);%设定旋转因子nk=n*k;%产生NN 维的nk整数矩阵WNnk=WN.nk;%用元素群运算求出W矩阵Xk=xn*WNnk;%求出离散傅里叶级数系数前两条语句无须解释。第三条语句把n转为列向量n，再与行向量k做矩阵乘，它产生的是一个整数矩阵：,实现频谱计算的程序(2),这个整数方阵恰好是算式中W矩阵的指数部分。第四条语句就产生了W矩阵，而最后一条语句就完成了DFT的全部运算。所以如果写成MATLAB子程序，它可以命名为DFT。在这5行DFT子程序前面，加上输入语句：xn=input(xn=);N=length(xn);在子程序后面，加上绘图语句：subplot(2,1,1),stem(xn,.);subplot(2,1,2),stem(abs(Xk),.)就得到本题的程序ag1009。运行ag1009，并按提示输入xn，即可在子图1上得到xn序列的波形，并在子图2上得到此信号的幅频特性。,例10.9 求频谱的数字例,例如输入序列为xn=ones(1,64),zeros(1,192)，则得到的时域波形及其频谱图如下。,不难看出，计算每一个DFT分量X(i)，需要做N次复数乘法和N-1次复数加法。而要完成整个DFT，求出X，则需要做NN次复数乘法和N(N-1)次复数加法。另外它占的数据存储量为NN。在MATLAB中每个数据都是以复数双精度格式保存的，所以总共要16NN 字节。,矩阵对科学计算的贡献,因此，当N较大时，MATLAB会拒绝执行这个程序，并警告内存不足。为此，在执行ag1009时，建议输入的xn长度不要超过256。此外，人们研究了很多种加快计算并节省内存的方法，统称为快速傅里叶变换(FFT)。从这里也可以看出，矩阵给我们提供了一个组织海量数据进行运算的最好方法，包括对上百万数据进行赋值、运算和绘图，只用了几行程序。比如其中给100多万个元素的DFT矩阵赋值，用的是1024元单列矩阵乘以单行矩阵的方法，这真是矩阵的奇妙点之一。线性代数的重要性之所以在近20年来可与微积分相妣美，这是一个重要原因。如果学线性代数而不涉及工程计算实践，那就无法真正体会线性代数的精髓所在。,10.10 最优FIR滤波器设计,设给定滤波器的预期幅频特性D()在i(i1,2,K)的K个频点上的值，若实际滤波器的脉冲响应的长度为N=2L+1，则其幅特性与预期特性相拟合的方程为：这K个方程联立。由于是L1个谐波的线性组合，其中待定参数d(n)有L1个。若KL1，为欠定方程，有无数解，无工程意义；若KL1，则方程组是适定的，恰好可以解出这些系数；若KL1，即方程数比未知数的数目多，成为超定方程组。那就只能找到近似地满足这些方程的最小二乘解d(n)，工程中常遇到的是这种问题。,超定条件下的方程,所以只考虑KL1的情况。eDPd其中：e为单列K元误差向量，D为预期幅特性在样本点列上的K元单列向量，d为L1元待定系数单列向量，而P则是由cos（ni）组成的K(L1)的系数矩阵。前面已经得到求d的最小二乘解的公式：,例10.10 FIR滤波器设计数字例,要求设计长度为N9(有5个待定系数)的FIR滤波器，使它在0之间的八个频点上逼近预期的低通幅特性D：0,0.33,0.67,1,1.5,2,2.5,3.14;D1,1,1,1,0,0,0,0;解：列出方程组的完整形式：,解此题目的程序ag1010,看来系数矩阵P的赋值比较麻烦，其实，回想求离散傅里叶变换时的做法，可以利用列矩阵ww1;w2;w8乘行矩阵0:4形成一个85矩阵，然后求余弦得到P。这可以用MATLAB语句Pcos(w*0:4)轻而易举地列出。因此用下列几行MATLAB程序ag1010就可以方便地完成最小二乘最优滤波器的设计。N=9;L=floor(N-1)/2);%列出待求系数序数%列出频率向量w=0,0.33,0.67,1,1.5,2,2.5,3.14;D=1,1,1,1,0,0,0,0;%列出预期幅特性向量P=cos(w*0:L);%列出P矩阵d=inv(P*P)*P*D%求出待求系数,程序ag1010运行的结果为,d=0.3981 0.6039 0.2137-0.1205-0.1506知道d(n)以后，就可以计算滤波器的频率特性进行校核。得到的幅频特性见图。因为L4，最高的谐波次数为4，在0之间幅特性摆动最多两个周期，达到峰值的次数应为四次，这从图上也看得很清楚。,10.11 信号流图的矩阵建模求解,任何复杂的结构的线性系统，其内部的信号传输关系，都可画成信号流图，在信号与系统、数字信号处理以及自动控制课程中所介绍的唯一解法是梅森公式。它没有证明，只能叫学生死记硬背，计算起来又非常繁琐，阶次高一点就极易出错，更无法用计算机辅助求解。所以在实际中极少有用处。利用MATLAB在机算矩阵问题上的巨大优势，用矩阵来给信号流图建模，可以得出简明的公式。设信号流图中有I个输入节点，N个中间和输出节点，它们分别代表输入信号ui（i=1,2,I）和系统状态xj（j=1,2,N）。信号流图代表它们之间的联结关系。用拉普拉斯算子表示后，任意状态xj可以表为ui和xj的线性组合：,信号流图的矩阵表示式,写成矩阵形式后为：其中：X为K维状态列向量,U为I维输入列向量,Q为NN维的联接矩阵，P为NI维的输入矩阵，Q和P的元素qji和pjk是各信号节点之间的传递系数。,信号流图矩阵的解,移项后，上式可写为由此得到各输入对各信号的传递函数公式：因此，系统的传递函数矩阵为W=(I-Q)-1 P，这个简明的公式就等价于梅森公式。只要写出P和Q，任何复杂系统的传递函数都可用这个简单的式子求出。存在的困难在于，连续系统的传递关系是用常数及拉普拉斯算子s组成的，离散系统的传递关系是用常数及z变换算子z或z-1组成的，但这个矩阵公式中用到的是普通的矩阵乘法和加法，无法应用于算子变量。用MATLAB的符号运算(Synbolic)工具箱解决了这个问题。从下面的数字例中可看出其做法。,例10.11 信号流图解法举例,如图，求以u为输入，x8为输出的传递函数。解：先列出信号流图方程组，把每一个信号写在等式左端，右面是其他信号的线性组合。信号间的传递关系先笼统地用符号G1,G2,来表示，得下列方程组 x1=ux2=x1-x3-x5x3=G1*x2x4=x3+x1-x5x5=G2*x4x6=x3+x5-x7x7=G3*x6x8=K*x7,例10.11的矩阵方程,用矩阵表示，可写成：以u为输入，X为输出的系统传递函数为W=(I-Q)-1 P它有8行，最后一行是W(8)=x8/u，即本题要求的结果。,程序ag1011,syms G1 G2 G3 K%定义字符自变量Q(3,2)=G1;%Q的第一条赋值语句右端必须是字符变量Q(2,1)=1;Q(2,3)=-1;Q(2,5)=-1;%列出连接矩阵Q(4,3)=1;Q(4,1)=1;Q(4,5)=-1;Q(5,4)=G2;Q(6,3)=1;Q(6,5)=1;Q(6,7)=-1;Q(7,6)=G3;Q(8,7)=K;Q(:,end+1)=zeros(max(size(Q),1)%加一个全零列P=1;0;0;0;0;0;0;0;I=eye(size(Q);W=(I-Q)P%求出完整的传递矩阵W8=W(8)%x8为输出的传递函数为其第八项W(8)pretty(W8)%给出便于阅读的形式,程序运行结果,W=1 1/(2*G1*G2+1+G2+G1)G1/(2*G1*G2+1+G2+G1)(2*G1+1)/(2*G1*G2+1+G2+G1)G2*(2*G1+1)/(2*G1*G2+1+G2+G1)(2*G1*G2+G2+G1)/(2*G1*G2+2*G1*G2*G3+1+G3+G2+G2*G3+G1+G1*G3)G3*(2*G1*G2+G2+G1)/(2*G1*G2+2*G1*G2*G3+1+G3+G2+G2*G3+G1+G1*G3)G3*K*(2*G1*G2+G2+G1)/(2*G1*G2+2*G1*G2*G3+1+G3+G2+G2*G3+G1+G1*G3)W8=G3*K*(2*G2*G1+G1+G2)/(2*G2*G1+2*G2*G1*G3+G3+G1+1+G2+G2*G3+G1*G3)写成美观形式为：读者如果熟悉梅森公式，可以用它算的结果进行比较，应该是完全相同的。同时，可以测定一下所化费的时间，也作一比较。,Gi是传递函数的情况,如果G1,G2,G3是自变量s的函数，那自变量中就不该有Gi出现，第一条语句应改为syms K s，再把其中的三条赋值语句改为：Q(3,2)=s/(s+1);Q(5,4)=3/(s+4);Q(7,6)=1/(s2+5*s+6);修改后的程序命名为ag1011a，运行此程序的结果为：这个例子是对连续系统的，用的算子变量是s，对离散系统，只要把算子变量换成时移算子z，或q=z-1，其他一切都同样处理即可，,10.12 自动控制系统的矩阵建模,众所周知，系统结构图与信号流图是等价的，所以上述信号流图的矩阵建模也适用于线性自动控制系统。在自控中，各条支路的运算函数一般比较复杂，因为这些支路往往代表一个部件，本身就是一个高阶的子系统。为了表示部件动态特性的方便，很少用信号流图，多用结构图。但这不要紧，矩阵建模的关键不在用什么图，而在于列写联立方程。用结构图同样可以列出系统的方程组。若系统比较简单，通常就可以用并联、串联和反馈三种简单情况的公式来求合成的传递函数，即使如此，笔算仍是令人头痛的事。若系统的节点数（也就是信号数）比较多，结构又比较复杂，按现在书上教的手工推导方法，往往就要把节点作移位，使系统变换成上述三种简单情况的叠加，这很容易出错，而且往往会因此改变一些信号的性质，失去其结构图上原有的重要物理意义。我们发现，用线性方程组的矩阵解法处理这个难题是最好的选择。,结构图的矩阵建模,首先根据系统结构图，自由选择若干个最关心的信号点组成向量，分别列出它们的方程，只要方程的右端出现的信号都在所选的信号向量X之中。这样就可写出方程组的矩阵形式，然后就用公式 求出对各点的传递函数。这样的不但解法可以保持各信号节点的原来意义，而且只要方程列写正确，推导过程的正确性将由计算软件来保证。,例10.12,对图所示的随动系统，求以驱动信号u1及测量干扰u2为输入，图中x1,x2,x3,x4四个信号节点为输出的传递函数，要求以美观方式显示系统的跟踪误差x1的表达式，说明它如何受u1,u2影响。解：列写方程，有：,本题的矩阵方程,写成矩阵形式为：左右对照，可得知系数矩阵Q和P的内容，根据就可以求出以u1,u2为输入，x1x4为输出的系统传递函数。这里的Q是44矩阵，当然要靠机算求(I-Q)的逆阵。x1点的信号可用求出。,程序ag1012,因此可列出解本题的程序ag1012如下：syms W1 W2 W3 W4 u1 u2 Q(1,4)=-W3;Q(2,1)=W1;%给Q的第一个元素赋值的是符号变量Q(3,2)=W2;Q(4,3)=1;P(2,1)=W4;P(1,1)=1;P(4,2)=1;%给P第一个赋值元素是符号变量W=inv(eye(4)-Q)*P%信号流图方程解x1=W(1,:)*u1;u2%计算各信号点在u1,u2作用下的输出x1pretty(W(1,1),pretty(W(1,2),%美观显示对两个传递函数pretty(x1)%美观显示输出误差项,程序运行结果,W=1/(1+w1*w2*w3)-w2*w3/(1+w1*w2*w3)*w4,-w3/(1+w1*w2*w3)w1/(1+w1*w2*w3)+1/(1+w1*w2*w3)*w4,-w1*w3/(1+w1*w2*w3)w2*w1/(1+w1*w2*w3)+w2/(1+w1*w2*w3)*w4,-w1*w2*w3/(1+w1*w2*w3)w2*w1/(1+w1*w2*w3)+w2/(1+w1*w2*w3)*w4,1/(1+w1*w2*w3)而 x1=(w2*w1/(1+w1*w2*w3)+w2/(1+w1*w2*w3)*w4)*u1-w1*w2*w3/(1+w1*w2*w3)*u2W(i,j)分别表示第j个输入对第i个输出的传递函数，W是42矩阵，说明这种方法的计算效率非常高，一次把8个传递函数都求出来了。不过其书写的形式太难读了，pretty是符号运算工具箱中的一个命令，它可以使表达式比较美观也易于阅读。得到：,程序ag1012a,用矩阵解法还可以对系统详细特性进行推导。例如，设此时可设拉普拉斯算子s为符号变量，W1,W2,W3为s的函数，放在程序ag1012的前部，把程序的前两句换成为下面5条语句，程序取名为ag1012a：clear,syms s u1 u2W1=1/(0.05*s+1);W2=2/(0.01*s2+0.1*s+1)/s;W3=1;W4=0.5*s/(0.001*s+1);Q(2,1)=W1;Q(1,4)=-W3;%先赋值符号变量W1.(以下同程序ag1012),程序ag1012a运行结果,运行此程序，将得出一个比较冗长的结果，因为Symbolic工具箱还不能进行近似处理，相近的因子在相减或相除中消不掉，人工近似处理后的结果为：手工推导这个结果恐费时很多，且出错概率很大，这是教育现代化（这里是用矩阵建模并用软件工具解题）给我们带来的好处，要充分利用。,