1、1第 2 章 符号计算所谓符号计算是指：解算数学表达式、方程不是在离散化的数值点上进行，而是凭借一系列恒等式，数学定理，通过推理和演绎，力求获得解析结果。这种计算建立在数值完全准确表达和推演严格解析的基础之上，因此所得结果是完全准确的。本书之所以把符号计算内容放在第 2 章，是出于以下考虑：一，相对于 MATLAB 的数值计算“引擎”和“函数库”而言，符号计算的“引擎”和“函数库”是独立的。二，在相当一些场合，符号计算解算问题的指令和过程，显得比数值计算更自然、更简明。三，大多数理工科的本科学生在学过高等数学和其他专业基础课以后，比较习惯符号计算的解题理念和模式。在编写本章时，作者在充分考虑符
2、号计算独立性的同时，还考虑了章节的自完整性。为此，本章不但全面地阐述符号计算，而且在最后一节还详细叙述了符号计算结果的可视化。 这样的安排，将使读者在阅读完本章后，就有可能运用 MATLAB 的符号计算能力去解决相当一些具体问题。2.1 符号对象和符号表达式2.1.1 符号对象的创建和衍生一 生成符号对象的基本规则二 符号数字【例 2.1-1】符号（类）数字与数值（类）数字之间的差异。a=pi+sqrt(5)sa=sym(pi+sqrt(5)Ca=class(a)Csa=class(sa)vpa(sa-a) a =5.3777sa =pi+sqrt(5)Ca =doubleCsa =syman
3、s =.138223758410852e-16 三 符号参数四 符号变量【例 2.1-2】用符号计算研究方程 的解。02wvzu（1）syms u v w z2Eq=u*z2+v*z+w;result_1=solve(Eq) %findsym(Eq,1) result_1 =-u*z2-v*zans =w （2）result_2=solve(Eq,z) result_2 =1/2/u*(-v+(v2-4*u*w)(1/2)1/2/u*(-v-(v2-4*u*w)(1/2) 【例 2.1-3】对独立自由符号变量的自动辨认。（1）syms a b x X Yk=sym(3);z=sym(c*sqr
4、t(delta)+y*sin(theta);EXPR=a*z*X+(b*x2+k)*Y; （2）findsym(EXPR)ans =X, Y, a, b, c, delta, theta, x, y （3）findsym(EXPR,1) ans =x （4）findsym(EXPR,2),findsym(EXPR,3) ans =x,yans =x,y,theta 【例 2.1-4】findsym 确定自由变量是对整个矩阵进行的。syms a b t u v x yA=a+b*x,sin(t)+u;x*exp(-t),log(y)+vfindsym(A,1) A = a+b*x, sin(t)
5、+u x*exp(-t), log(y)+vans =x 32.1.2 符号计算中的算符2.1.3 符号计算中的函数指令2.1.4 符号对象的识别【例 2.1-5】数据对象及其识别指令的使用。（1）cleara=1;b=2;c=3;d=4;Mn=a,b;c,dMc=a,b;c,dMs=sym(Mc) Mn =1 23 4Mc =a,b;c,dMs = a, b c, d （2）SizeMn=size(Mn)SizeMc=size(Mc)SizeMs=size(Ms) SizeMn =2 2SizeMc =1 9SizeMs =2 2 （3）CMn=class(Mn)CMc=class(Mc)C
6、Ms=class(Ms) CMn =doubleCMc =charCMs =sym （4）isa(Mn,double)isa(Mc,char)isa(Ms,sym) ans =1ans =1ans =1 （5）whos Mn Mc Ms Name Size Bytes Class Attributes4Mc 1x9 18 char Mn 2x2 32 double Ms 2x2 312 sym 2.2 符号数字及表达式的操作2.2.1 数值数字与符号数字之间的转换一 数值数字向符号数字的转换二 符号数字向双精度数字转换2.2.2 符号数字的任意精度计算【例 2.2-1】digits, vpa
7、指令的使用。digitsp0=sym(1+sqrt(5)/2)pr=sym(1+sqrt(5)/2)%pd=sym(1+sqrt(5)/2,d)%e32r=vpa(abs(p0-pr) )e16=vpa(abs(p0-pd),16) e32d=vpa(abs(p0-pd)Digits = 32p0 =(1+sqrt(5)/2pr =7286977268806824*2(-52)pd =1.6180339887498949025257388711907e32r =.543211520368251e-16e16 =0.e32d =.543211520368251e-16 2.2.3 符号表达式的基
8、本操作【例 2.2-2】简化 。32816xfsyms xf=(1/x3+6/x2+12/x+8)(1/3);g1=simple(f)g2=simple(g1) g1 =(2*x+1)/xg2 =2+1/x 52.2.4 表达式中的置换操作一 子表达式置换操作【例 2.2-3】对符号矩阵 进行特征向量分解。dcbaclear allsyms a b c d WV,D=eig(a b;c d)RVD,W=subexpr(V;D,W) V = -(1/2*d-1/2*a-1/2*(d2-2*a*d+a2+4*b*c)(1/2)/c, -(1/2*d-1/2*a+1/2*(d2-2*a*d+a2+4
9、*b*c)(1/2)/c 1, 1D = 1/2*d+1/2*a+1/2*(d2-2*a*d+a2+4*b*c)(1/2), 0 0, 1/2*d+1/2*a-1/2*(d2-2*a*d+a2+4*b*c)(1/2)RVD = -(1/2*d-1/2*a-1/2*W)/c, -(1/2*d-1/2*a+1/2*W)/c 1, 1 1/2*d+1/2*a+1/2*W, 0 0, 1/2*d+1/2*a-1/2*WW =(d2-2*a*d+a2+4*b*c)(1/2) 二 通用置换指令【例 2.2-4】用简单算例演示 subs 的置换规则。（1）syms a x;f=a*sin(x)+5 f =a
10、*sin(x)+5 （2）f1=subs(f,sin(x),sym(y) %class(f1) f1 =a*y+5ans =sym （3）f2=subs(f,a,x,2,sym(pi/3) %class(f2) f2 =3(1/2)+5ans =sym 6（4）f3=subs(f,a,x,2,pi/3) %class(f3) f3 =6.7321ans =double （5）f4=subs(subs(f,a,2),x,0:pi/6:pi) %class(f4) f4 =5.0000 6.0000 6.7321 7.0000 6.7321 6.0000 5.0000ans =double （6）
11、f5=subs(f,a,x,0:6,0:pi/6:pi) %class(f5) f5 =5.0000 5.5000 6.7321 8.0000 8.4641 7.5000 5.0000ans =double 2.3 符号微积分2.3.1 极限和导数的符号计算【例 2.3-1】试求 。kxx1limsyms x kLim_f=limit(1-1/x)(k*x),x,inf) Lim_f =exp(-k) 【例 2.3-2】求 求 ， ， 。xttaflncos3df2tfdxfsyms a t x;f=a,t3;t*cos(x), log(x);df=diff(f)dfdt2=diff(f,t,
12、2)dfdxdt=diff(diff(f,x),t) df = 0, 0 -t*sin(x), 1/xdfdt2 = 0, 6*t 0, 0dfdxdt = 0, 0 -sin(x), 0 7【例 2.3-3】求 的 Jacobian 矩阵 。)sin(co),(21212xexxf 23121xfxfsyms x1 x2;f=x1*exp(x2);x2;cos(x1)*sin(x2);v=x1 x2;fjac=jacobian(f,v) fjac = exp(x2), x1*exp(x2) 0, 1 -sin(x1)*sin(x2), cos(x1)*cos(x2) 【例 2.3-4】 ，求
13、 ， 。xfsin)()0(xffx（1）clearsyms xsyms d positivef_p=sin(x); %df_p=limit(subs(f_p,x,x+d)-f_p)/d,d,0) % df_p0=limit(subs(f_p,x,d)-subs(f_p,x,0)/d,d,0) % df_p =cos(x)df_p0 =1 （2）f_n=sin(-x);df_n=limit(f_n-subs(f_n,x,x-d)/d,d,0) % df_n0=limit(subs(f_n,x,0)-subs(f_n,x,-d)/d,d,0) % df_n =-cos(x)df_n0 =-1 （
14、3）f=sin(abs(x);dfdx=diff(f,x) % dfdx =cos(abs(x)*abs(1,x) （4）xn=-3/2*pi:pi/50:0;xp=0:pi/50:3/2*pi;xnp=xn,xp(2:end);hold onplot(xnp,subs(f,x,xnp),k,LineWidth,2.5) % plot(xn,subs(df_n,x,xn),-r,LineWidth,2.5)plot(xp,subs(df_p,x,xp),:r,LineWidth,2.5)legend(char(f),char(df_n),char(df_p),Location,NorthEas
15、t)% grid onxlabel(x)hold off 8图 2.3-1 函数及其导函数【例 2.3-5】设 ，求 。yxsin)icos(dx（1）clearsyms xg=sym(cos(x+sin(y(x)=sin(y(x) % dgdx=diff(g,x) % g =cos(x+sin(y(x)=sin(y(x)dgdx =-sin(x+sin(y(x)*(1+cos(y(x)*diff(y(x),x) = cos(y(x)*diff(y(x),x) （2）disp(maple(isolate,dgdx,diff(y(x),x) % diff(y(x),x) = sin(x+sin(
16、y(x)/(-sin(x+sin(y(x)*cos(y(x)-cos(y(x) 【例 2.3-6】求 在 处展开的 8 阶 Maclaurin 级数。xef)(0syms xr=taylor(x*exp(x),9,x,0) %忽略 9 阶及 9 阶以上小量的展开 r =x+x2+1/2*x3+1/6*x4+1/24*x5+1/120*x6+1/720*x7+1/5040*x8 【例 2.3-7】求 在 处的截断 8 阶小量的 Taylor 展开近似。)sin(2yx0,TL1=maple(mtaylor(sin(x2+y),x=0,y=0,8)class(TL1) TL1 =y+x2-1/6*
17、y3-1/2*x2*y2+1/120*y5-1/2*x4*y+1/24*x2*y4-1/6*x6-1/5040*y7+1/12*x4*y3ans =char 92.3.2 序列/级数的符号求和【例 2.3-8】求 ， 。103,ttk12)1(,(kksyms k t;f1=t k3;f2=1/(2*k-1)2,(-1)k/k;s1=simple(symsum(f1)s2=simple(symsum(f2,1,inf) s1 = 1/2*t*(t-1), k3*ts2 = 1/8*pi2, -log(2) 2.3.3 符号积分【例 2.3-9】求 。dx1clearsyms xf=sqrt(1
18、+x)/x)/xs=int(f,x)s=simple(simple(s) f =(1+x)/x)(1/2)/xs =(1+x)/x)(1/2)/x*(-2*(x2+x)(3/2)+2*(x2+x)(1/2)*x2+log(1/2+x+(x2+x)(1/2)*x2)/(1+x)*x)(1/2)s =log(1/2+x+(1+x)*x)(1/2)-2*(1+x)*x)(1/2)/x 【例 2.3-10】求 。dxxbasin12syms a b x;f=a*x,b*x2;1/x,sin(x);disp(The integral of f is);pretty(int(f) The integral
19、 of f is 2 31/2 a x 1/3 b x log(x) -cos(x) 【例 2.3-11】求积分 。2 1 22)(xydzyxsyms x y zF2=int(int(int(x2+y2+z2,z,sqrt(x*y),x2*y),y,sqrt(x),x2),x,1,2)VF2=vpa(F2) F2 =1610027357/6563700+64/225*2(3/4)-6072064/348075*2(1/2)+14912/4641*2(1/4)VF2 =224.92153573331143159790710032805 10【例 2.3-12】求阿基米德（Archimedes）
20、螺线 在 到 间的曲线长)0(,ar度函数，并求出 时的曲线长度。2,1a（1）syms a r theta phi positivex=r*cos(theta);x=subs(x,r,a*theta);y=r*sin(theta);y=subs(y,r,a*theta);dLdth=sqrt(diff(x,theta)2+diff(y,theta)2);L=simple(int(dLdth,theta,0,phi) L =1/2*a*(phi*(phi2+1)(1/2)-log(-phi+(phi2+1)(1/2) （2）L_2pi=subs(L,a,phi,sym(1,2*pi)L_2pi
21、_vpa=vpa(L_2pi) L_2pi =pi*(1+4*pi2)(1/2)+1/2*asinh(2*pi)L_2pi_vpa =21.256294148209098800702512272566 （3）L1=subs(L,a,sym(1);ezplot(L1*cos(phi),L1*sin(phi),0,2*pi)grid onhold onx1=subs(x,a,sym(1);y1=subs(y,a,sym(1);h1=ezplot(x1,y1,0,2*pi);set(h1,Color,r,LineWidth,5)title( )legend(螺线长度 -幅角曲线, 阿基米德螺线)ho
22、ld off 图 2.3-2 阿基米德螺线（粗红）和螺线长度函数（细蓝）2.4 微分方程的符号解法112.4.1 符号解法和数值解法的互补作用2.4.2 求微分方程符号解的一般指令2.4.3 微分方程符号解示例【例 2.4-1】求 的解。dxtytx,S=dsolve(Dx=y,Dy=-x)disp(blanks(12),x,blanks(21),y),disp(S.x,S.y) S = x: 1x1 symy: 1x1 symx y -C1*cos(t)+C2*sin(t), C1*sin(t)+C2*cos(t) 【例 2.4-2】图示微分方程 的通解和奇解的关系。2)(yxy=dsolv
23、e(Dy)2-x*Dy+y=0,x) % clf,hold on,ezplot(y(1),-6,6,-4,8,1) %cc=get(gca,Children); %set(cc,Color,r,LineWidth,5) %for k=-2:0.5:2;ezplot(subs(y(2),C1,k),-6,6,-4,8,1);end, %hold off,title(fontname隶书fontsize16通解和奇解) y =1/4*x2-C12+x*C1 图 2.4-1 通解和奇解曲线【例 2.4-3】求解两点边值问题： 。0)5(,)1,32yxy（1）y=dsolve(x*D2y-3*Dy=
24、x2,y(1)=0,y(5)=0,x) y =31/468*x4-1/3*x3+125/468 （2）12ezplot(y,-1,6)hold onplot(1,5,0,0,.r,MarkerSize,20)text(1,1,y(1)=0)text(4,1,y(5)=0)title(x*D2y-3*Dy=x2, y(1)=0,y(5)=0)hold off 图 2.4-2 两点边值问题的解曲线2.5 符号变换和符号卷积2.5.1 Fourier 变换及其反变换【例 2.5-1】求 的 Fourier 变换。01)(ttf（1）syms t wut=heaviside(t);UT=fourier
25、(ut) UT =pi*dirac(w)-i/w （2）Ut=ifourier(UT,w,t) Ut =heaviside(t) 【例 2.5-2】根据 Fourier 变换定义，用积分指令求方波脉冲的 Fourier 变换。elstAy2/0（1）syms A t wsyms tao positive13yt=heaviside(t+tao/2)-heaviside(t-tao/2);Yw=fourier(A*yt,t,w) Yw =2*A/w*sin(1/2*tao*w) （2）Yt=ifourier(Yw,w,t) Yt =A*(heaviside(t+1/2*tao)-heavisid
26、e(t-1/2*tao) （3）yt3=subs(yt,tao,3)Yw3=subs(Yw,A,tao,1,3)subplot(2,1,1)Ht=ezplot(yt3,-3,3);set(Ht,Color,r,LineWidth,3)subplot(2,1,2),ezplot(Yw3) yt3 =heaviside(t+3/2)-heaviside(t-3/2)Yw3 =2/w*sin(3/2*w) -3 -2 -1 0 1 2 300.51theaviside(t+3/2)-heaviside(t-3/2)-6 -4 -2 0 2 4 6-10123w2/w sin(3/2 w)图 2.5-
27、1 时域方波及其 Fourier 变换【例 2.5-3】求 的 Fourier 变换，在此 是参数， 是时间变量。xtetfxt0)()( xtsyms t x w;ft=exp(-(t-x)*heaviside(t-x);F1=simple(fourier(ft,t,w)F2=simple(fourier(ft)F3=simple(fourier(ft,t) F1 =1/(1+i*w)/exp(i*x*w)F2 =i*exp(-i*t*w)/(i+w)F3 =-exp(-t*(2+i*t)/(-1+i*t) 142.5.2 Laplace 变换及其反变换【例 2.5-4】求 的 Laplac
28、e 变换。tatebeu2sin)()(syms t s;syms a b positive % Dt=dirac(t-a);Ut=heaviside(t-b);Mt=Dt,Ut;exp(-a*t)*sin(b*t),t2*exp(-t);MS=laplace(Mt,t,s) MS = exp(-s*a), exp(-s*b)/s 1/b/(s+a)2/b2+1), 2/(s+1)3 【例 2.5-5】验证 Laplace 时移性质： 。0)()(000tfLetUtfLstsyms t ssyms t0 positiveft=sym(f(t-t0)*heaviside(t-t0)FS=lap
29、lace(ft,t,s)FS_t=ilaplace(FS,s,t) ft =f(t-t0)*heaviside(t-t0)FS =exp(-s*t0)*laplace(f(t),t,s)FS_t =f(t-t0)*heaviside(t-t0) 2.5.3 Z 变换及其反变换【例 2.5-6】求序列 的 Z 变换，并用反变换验算。fnn().)026150（1）syms nDelta=sym(charfcn0(n); % 定义离散单位函数 01)(nD0=maple(char(subs(Delta,n,0) %计算 0D15=maple(char(subs(Delta,n,15) %计算 )5
30、( D0 =1D15 =0 （2）syms zfn=2*Delta+6*(1-(1/2)n)FZ=simple(ztrans(fn,n,z);disp(FZ = )pretty(FZ)FZ_n=iztrans(FZ,z,n) fn =2*charfcn0(n)+6-6*(1/2)nFZ = 1524 z + 2-22 z - 3 z + 1FZ_n =2*charfcn0(n)+6-6*(1/2)n 2.5.4 符号卷积【例 2.5-7】已知系统冲激响应 ，求 输入下的输出响htTeUtt()()/1uteUt()()应。syms T t taout=exp(-t);ht=exp(-t/T)/
31、T;uh_tao=subs(ut,t,tao)*subs(ht,t,t-tao);yt=simple(simple(int(uh_tao,tao,0,t) yt =-(exp(-t)-exp(-t/T)/(T-1) 【例 2.5-8】采用 Laplace 变换和反变换求上例的输出响应。syms syt=ilaplace(laplace(ut,t,s)*laplace(ht,t,s),s,t);yt=simple(yt) yt =(-exp(-t)+exp(-t/T)/(T-1) 【例 2.5-9】求函数 和 的卷积。utUt()()1hteUt()syms taot=sym(t,positiv
32、e);ut=heaviside(t)-heaviside(t-1);ht=t*exp(-t);yt=int(subs(ut,t,tao)*subs(ht,t,t-tao),tao,0,t)yt=collect(yt,heaviside(t-1) yt =1+t*heaviside(t-1)*exp(1-t)-heaviside(t-1)+(-t-1)/exp(t)yt =(exp(1-t)*t-1)*heaviside(t-1)+1+(-t-1)/exp(t) 2.6 符号矩阵分析和代数方程解2.6.1 符号矩阵分析【例 2.6-1】求矩阵 的行列式、逆和特征根。Aa12syms a11 a1
33、2 a21 a22A=a11,a12;a21,a22DA=det(A)IA=inv(A)16EA=eig(A) A = a11, a12 a21, a22DA =a11*a22-a12*a21IA = a22/(a11*a22-a12*a21), -a12/(a11*a22-a12*a21) -a21/(a11*a22-a12*a21), a11/(a11*a22-a12*a21)EA =1/2*a11+1/2*a22+1/2*(a112-2*a11*a22+a222+4*a12*a21)(1/2)1/2*a11+1/2*a22-1/2*(a112-2*a11*a22+a222+4*a12*a
34、21)(1/2) 【例 2.6-2】著名的 Givens 旋转（变换） 对矩阵ttGcosini的旋转作用。231A（1）syms tA=sym(sqrt(3)/2,1/2;1/2,sqrt(3)/2)G=cos(t),-sin(t);sin(t),cos(t);GA=G*A A = sqrt(3/4), 1/2 1/2, sqrt(3/4)GA = 1/2*cos(t)*3(1/2)-1/2*sin(t), 1/2*cos(t)-1/2*sin(t)*3(1/2) 1/2*sin(t)*3(1/2)+1/2*cos(t), 1/2*sin(t)+1/2*cos(t)*3(1/2) （2）An
35、=subs(GA,t,110/180*pi);Op=0;0;Ad=double(A);v1=Op,Ad(:,1);v2=Op,Ad(:,2);u1=Op,An(:,1);u2=Op,An(:,2);plot(v1(:,1),v1(:,2),-k,v2(:,1),v2(:,2),b)axis(-1,1,-1,1),axis squarehold onhu=plot(u1(:,1),u1(:,2),-k,u2(:,1),u2(:,2),b);set(hu,LineWidth,4)title(Givens Rotation)legend( v1, v2, u1, u2,Location,South)
36、hold off grid on 17-1 -0.5 0 0.5 1-1-0.8-0.6-0.4-0.200.20.40.60.81Givens Rotationv1 v2 u1u2图 2.6-1 Givens 旋转的几何意义【例 2.6-3】求 线性dnpqdpqdnpqnd21481,方程组的解。（1）A=sym(1 1/2 1/2 -1;1 1 -1 1;1 -1/4 -1 1;-8 -1 1 1);b=sym(0;10;0;1)X1=Ab b =01001X1 =1889 （2）S=solve(d+n/2+p/2-q,d+n-p+q-10,d-n/4-p+q,-8*d-n+p+q-1)
37、;disp( d, n, p, q)disp(S.d,S.n,S.p,S.q) d n p q 1, 8, 8, 9 2.6.2 线性方程组的符号解2.6.3 一般代数方程组的解【例 2.6-4】求方程组 ， 关于 的解。uyvzw20yz0zy,S=solve(u*y2+v*z+w=0,y+z+w=0,y,z)disp(S.y),disp(S.y),disp(S.z),disp(S.z) S = y: 2x1 symz: 2x1 sym18S.y-1/2/u*(-2*u*w-v+(4*u*w*v+v2-4*u*w)(1/2)-w-1/2/u*(-2*u*w-v-(4*u*w*v+v2-4*u
38、*w)(1/2)-wS.z1/2/u*(-2*u*w-v+(4*u*w*v+v2-4*u*w)(1/2)1/2/u*(-2*u*w-v-(4*u*w*v+v2-4*u*w)(1/2) 【例 2.6-5】solve 指令求 ， ， 构成的“欠qpnd210pdpndq4定”方程组解。syms d n p qeq1=d+n/2+p/2-q;eq2=n+d+q-p-10;eq3=q+d-n/4-p;S=solve(eq1,eq2,eq3,d,n,p,q);disp( S.d, S.n, S.p, S.q)disp(S.d,S.n,S.p,S.q) Warning: 3 equations in 4
39、variables. In solve at 113In sym.solve at 49S.d S.n S.p S.q d, 8, 4*d+4, 3*d+6 【例 2.6-6】求 的解。2)(xclear all,syms x;s=solve(x+2)x=2,x) s =.69829942170241042826920133106081 2.7 符号传递函数求取的状态方程法2.7.1 结构框图的代数状态方程解法【例 2.7-1】求图 2.7-1 所示某三环系统的传递函数。本例演示：（ A）系统“代数状态方程”的建立；（B）根据代数状态方程求系统的传递函数；（C ）编写 M 码时，系统矩阵的输入
40、采用“全元素赋值法”；（D ）sort 指令对符号表达式的排序功能。图 2.7-1 三环系统的结构框图（1）在结构框图上标识状态变量（3）据代数状态方程求传递函数的理论演绎19（4）代数状态方程法计算传递函数的 M 码syms G1 G2 G3 G4 H1 H2 H3A= 0, 0, 0, 0, 0, 0, -G1;G2, 0, 0, 0, 0, -G2, 0;0, G3, 0, 0, G3, 0, 0;0, 0, G4, 0, 0, 0, 0;0, 0, 0, H2, 0, 0, 0;0, 0, 0, H1, 0, 0, 0;0, 0, 0, H3, 0, 0, 0;b= G1; 0; 0;
41、 0; 0; 0; 0;c= 0, 0, 0, 1, 0, 0, 0;Y2Ua=c*(eye(size(A)-A)b);NN,D0=numden(Y2Ua);D0DD=sort(D0)disp(blanks(5),传递函数 Y2Ua 为)pretty(NN/DD) D0 =G3*G2*G1*H3*G4+G3*G2*H1*G4+1-G3*H2*G4DD =G1*G2*G4*H3*G3+G2*G4*G3*H1-G4*H2*G3+1传递函数 Y2Ua 为G3 G4 G2 G1-G1 G2 H3 G4 G3 + G2 H1 G4 G3 - H2 G4 G3 + 1 2.7.2 信号流图的代数状态方程解
42、法【例 2.7-2】作为比较，画出图 2.7-1 所示 结构框图的等价信号流图，并据此信号流图运用“代数状态方程”求系统的传递函数。本例演示：（A ）信号流图的代数状态方程的建立；（B）根据代数状态方程求传递函数；（ C）在 ， ，10sG12s， ， ， ， 的情况下，求取4123sG61sG12sH3H参数具体化的传递函数。（1）根据图 2.7-1 画出相应的信号流图图 2.7-2 三环系统的信号流图（2）代数状态方程法求取信号流图传递函数的数学原理20（3）实现以上算法的 M 码syms G1 G2 G3 G4 H1 H2 H3A= 0, 0, 0, 0, -H3;G1, 0, 0, 0
43、, -H1;0, G2, 0, 0, H2;0, 0, G3, 0, 0;0, 0, 0, G4, 0;b= 1; 0; 0; 0; 0;c= 0, 0, 0, 0, 1;Y2Ub=c*(eye(size(A)-A)b);NN,DD=numden(Y2Ub);DD=sort(DD);disp(blanks(5),传递函数 Y2Ub 为)pretty(NN/DD) 传递函数 Y2Ub 为G3 G4 G2 G1-G1 G2 H3 G4 G3 + G2 H1 G4 G3 - H2 G4 G3 + 1 （4）方块参数具体化时的传递函数syms sY2Uc=simple(subs(Y2Ub,G1,G2,
44、G3,G4,H1,H2,H3,100/(s+10),1/(s+1),(s+1)/(s2+4*s+4),(s+1)/(s+6),(2*s+12)/(s+1),(s+1)/(s+2),1);NN,DD=numden(Y2Uc);NN=expand(NN);DD=sort(DD);disp(参数具体化的传递函数 Y2Uc 为)pretty(NN/DD) 参数具体化的传递函数 Y2Uc 为2100 s + 300 s + 200-5 4 3 2s + 21 s + 157 s + 663 s + 1301 s + 910 2.7.3 多输入多输出系统传递矩阵的求取【例 2.7-3】运用“代数状态方程法
45、 ”求图 2.7-3 所示“2 输入 2 输出”系统的传递矩阵。本例演示：（A）多输入多输出系统传递矩阵的 “代数状态方程”求取法；（B）在书写M 码时，系统状态矩阵创建的“关联元素赋值法”；（C）对于复杂符号表达式，如何采用公因式简化表达；（D）所求复杂系统传递函数的正确性检验。21图 2.7-3 二输入二输出复杂系统的结构框图（1）多输入多输出系统的数学描述（2）求取传递矩阵的状态方程法的 M 码实现syms G1 G2 G3 G4 G5 G6 H1 H2 H3syms F1 F2 F3 F4 F5 F6 Q1 Q2syms WA=sym(zeros(17,17);B=sym(zeros(17,2);C=sym(zeros(2,17);A(1,7)= -G1;B(1,1)=G1;A(2,1)
