《Matlab软件基础及应用》.doc

第6章 符号计算6.1 符号对象和符号表达式在MATLAB中，数值和数值变量用于数值的存贮和各种数值计算。而符号常数、符号变量、符号函数、符号操作等则是用来形成符号表达式，严格按照代数、微积分等课程中的规则、公式进行运算，并尽可能给出解析结果。6.1.1 符号对象的生成和使用f=sym(arg) 把数字、字符串或表达式arg转换成为符号对象ff=sym(argn,flagn) 把数字、字符串或表达式argn转换成为flagn格式的符号对象f flagn取：d最接近的十进制浮点精确表示； r最接近有理表示；argv=sym(argv,flagv) 按flagv指定的要求把字符串argv定义为符号对象argvsyms(argv1,argv2,argvk) 把字符argv1,argv2,argvk定义为基本符号对象syms argv1 argv2 argv3 上述格式的简单形式【例6.1-1】符号常数形成中的差异a1=1/3,pi/7,sqrt(5),pi+sqrt(5)a2=sym(1/3,pi/7,sqrt(5),pi+sqrt(5)a3=sym(1/3,pi/7,sqrt(5),pi+sqrt(5)a23=a2-a3 a1 = 0.3333 0.4488 2.2361 5.3777a2 = 1/3, pi/7, sqrt(5), 6054707603575008*2(-50)a3 = 1/3, pi/7, sqrt(5), pi+sqrt(5)a23 = 0, 0, 0, 189209612611719/35184372088832-pi-5(1/2) 【例6.1-2】把字符表达式转换为符号变量y=sym(2*sin(x)*cos(x)y=simple(y) y =2*sin(x)*cos(x)y =sin(2*x) 【例6.1-3】用符号计算验证三角等式。syms fai1 fai2y=simple(sin(fai1)*cos(fai2)-cos(fai1)*sin(fai2) y =sin(fai1-fai2) 【例6.1-4】求矩阵的行列式值、逆和特征根syms a11 a12 a21 a22A=a11,a12;a21,a22DA=det(A),IA=inv(A),EA=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*a21)(1/2) 【例6.1-5】验证积分。syms A t tao wyf=int(A*exp(-i*w*t),t,-tao/2,tao/2);Yf=simple(yf) Yf =2*A*sin(1/2*tao*w)/w 62 符号计算中的算符和基本函数（1） 基本运算符矩阵算符 * / 元素对元素算符 .* ./ . .621 识别对象类别的指令【例6.1-6】数据对象及其识别指令的使用。clear,a=1;b=2;c=3;d=4;Mn=a,b;c,dMc=a,b;c,dMs=sym(Mc) Mn = 1 2 3 4Mc =a,b;c,dMs = a, b c, d (1) 判断矩阵的大小SizeMn=size(Mn),SizeMc=size(Mc),SizeMs=size(Ms) SizeMn = 2 2SizeMc = 1 9SizeMs = 2 2 (2)用class获得每种矩阵的类别CMn=class(Mn),CMc=class(Mc),CMs=class(Ms) CMn =doubleCMc =charCMs =sym (3) 用isa判断矩阵的类别isa(Mn,double),isa(Mc,char),isa(Ms,sym) ans = 1ans = 1ans = 1 (4)利用whos观察内存变量的类别和其它属性whos Mn Mc Ms Name Size Bytes Class Mc 1x9 18 char array Mn 2x2 32 double array Ms 2x2 408 sym objectGrand total is 21 elements using 458 bytes 6.2.2 符号表达式中自由变量的确定findsym(EXPR) 确认表达式EXPR中所有自由符号变量。findsym(EXPR，N) 从表达式EXPR中确认出靠x最近的N个独立变量。【例6.1-7】对独立自由符号变量的自动辨认。syms a b x X Y;k=sym(3);z=sym(c*sqrt(delta)+y*sin(theta);EXPR=a*z*X+(b*x2+k)*Y; findsym(EXPR)ans =X, Y, a, b, c, delta, theta, x, y findsym(EXPR,1) ans =x findsym(EXPR,2),findsym(EXPR,3) ans =x,yans =x,y,theta 【例6.1-8】findsym确定自由变量是对整个矩阵进行的。syms a b t u v x y;A=a+b*x,sin(t)+u;x*exp(-t),log(y)+vfindsym(A,1) A = a+b*x, sin(t)+u x*exp(-t), log(y)+vans =x 63 符号对象的操作和转换631 符号表达式的操作 简化simple(expr)【例6.2-1】简化syms xf=(1/x3+6/x2+12/x+8)(1/3);g1=simple(f)g2=simple(g1) g1 =(2*x+1)/xg2 =2+1/x 合并同类项 collect（expr，v）因式分解 factor（expr）6.3.2 置换操作6.3.2.1 子表达式置换操作 rs ssub=subexpr(S,ssub) 运用符号变量ssub置换子表达式，重写S为rs【例6.2-2】把复杂表达式中所含的多个相同子表达式用一个符号代替，使表达简洁。clear allsyms a b c d WV,D=eig(a b;c d);RVD,W=subexpr(V;D,W) 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) 6.3.2.2 通用置换指令 RES=subs(ES,old,new) 用new置换ES中的old后产生RES RES=subs(ES,new) 用new置换ES中的自由变量后产生RES【例6.2-3】用简单算例演示subs的置换规则。syms a x;f=a*sin(x)+5; （1）符号变量置换f1=subs(f,sin(x),sym(y) f1 =a*y+5 （2）符号常数置换f2=subs(f,a,x,2,sym(pi/3) f2 =3(1/2)+5 （3）双精度数值置换f3=subs(f,a,x,2,pi/3) f3 = 6.7321 （4）数值数组置换一（取a=2,x=0:pi/6:pi）f4=subs(subs(f,a,2),x,0:pi/6:pi) f4 = 5.0000 6.0000 6.7321 7.0000 6.7321 6.0000 5.0000 （4）数值数组置换二(取a=0:6,x=0:pi/6:pi）f5=subs(f,a,x,0:6,0:pi/6:pi) f5 = 5.0000 5.5000 6.7321 8.0000 8.4641 7.5000 5.0000 6.3.2.3 符号数值精度控制和任意精度计算 【例6.2-4】指令使用演示。p0=sym(1+sqrt(5)/2);p1=sym(1+sqrt(5)/2)pd=double(p0)e01=vpa(abs(p0-p1)e0d=vpa(abs(p0-pd)p1 =7286977268806824*2(-52)pd = 1.6180e01 =.543211520368250e-16e0d =.543211520368250e-16 p2=vpa(p0)e02=vpa(abs(p0-p2),64) p2 = 1.6180339887498948482045868343657 e02 = .61882279690820194237137864551377e-31 digits Digits = 32 6.3.3 符号对象与其它数据对象间的转换6.4 符号微积分6.4.1 符号序列的求和 s=symsum(f,v,a,b) 求通式f在指定变量v取遍a,b中所有整数时的和【例6.3-1】求，syms 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) 642 符号微分和矩阵【例6.3-2】求、和syms a t x;f=a,t3;t*cos(x), log(x);df=diff(f)dfdt2=diff(f,t,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 【例6.3.2-2】求的矩阵。syms 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) 643 符号积分intf=int(f,v) 给出f对指定变量v的不定积分intf=int(f,v,a,b) 给出f对指定变量v的定积分【例6.3-4】求。syms a b x;f=a*x,b*x2;1/x,sin(x);disp(The integral of f is);pretty(int(f) The integral of f is 2 3 1/2 a x 1/3 b x log(x) -cos(x) 【例6.3-5】求。F1=int(1/log(t),t,0,x) F1 =-Ei(1,-log(x) x=0.5:0.1:0.9F115=-mfun(Ei,1,-log(x) x = 0.5000 0.6000 0.7000 0.8000 0.9000F115 = -0.3787 -0.5469 -0.7809 -1.1340 -1.7758 syms 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-6072064/348075*2(1/2)+14912/4641*2(1/4)+64/225*2(3/4)VF2 =224.92153573331143159790710032805 644 符号卷积【例5.3-7】本例演示卷积的时域积分法：已知系统冲激响应，求输入下的输出响应。syms T t positivesyms taout=exp(-t);ht=exp(-t/T)/T;uh_tao=subs(ut,t,tao)*subs(ht,t,t-tao);yt=int(uh_tao,tao,0,t);yt=simple(yt) yt =(-exp(-t)+exp(-t/T)/(T-1) 【例6.3-8】本例演示通过变换和反变换求取卷积。系统冲激响应、输入同上例，求输出。syms syt=ilaplace(laplace(ut,t,s)*laplace(ht,t,s),s,t);yt=simple(yt) yt =-(exp(-t)-exp(-t/T)/(T-1) 【例6.3-9】求函数和的卷积。syms tao;t=sym(t,positive);ut=sym(Heaviside(t)-Heaviside(t-1);ht=t*exp(-t);yt53=int(subs(ut,t,tao)*subs(ht,t,t-tao),tao,0,t);yt53=collect(yt53,Heaviside(t-1) yt53 =(-1+exp(1-t)*t)*Heaviside(t-1)+1+(-t-1)*exp(-t) 6.4.5 符号积分变换6.4.5.1 Fourier变换及其反变换【例 6.4-1】求的Fourier变换。syms t wut=sym(Heaviside(t);UT=fourier(ut)UTC=maple(convert,UT,piecewise,w)UTS=simple(UT) UT =pi*Dirac(w)-i/wUTC =PIECEWISE(undefined, w = 0,0, otherwise)UTS =-i/w Ut=ifourier(UT,w,t)Uts=ifourier(UTS,w,t)Ut =1/2+1/2*Heaviside(t)-1/2*Heaviside(-t)Uts =1/2*Heaviside(t)-1/2*Heaviside(-t) 【例6.4-2】用fourier指令求方波脉冲的Fourier变换。（1）syms A t wsyms tao positiveyt=sym(Heaviside(t+tao/2)-Heaviside(t-tao/2);Yw=fourier(A*yt,t,w)Ywc=maple(convert,Yw,piecewise,w)Yws=simple(Yw) Yw =A*(exp(1/2*i*tao*w)*(pi*Dirac(w)-i/w)-exp(-1/2*i*tao*w)*(pi*Dirac(w)-i/w)Ywc =-i*A*exp(1/2*i*tao*w)/w+i*A*exp(-1/2*i*tao*w)/wYws =2*A*sin(1/2*tao*w)/w （2）Yt=ifourier(Yw,w,t)Yst=ifourier(Yws,w,t) Yt =-1/2*A*(-Heaviside(t+1/2*tao)+Heaviside(-t-1/2*tao)+Heaviside(t-1/2*tao)-Heaviside(-t+1/2*tao)Yst =-1/2*A*(-Heaviside(t+1/2*tao)+Heaviside(-t-1/2*tao)+Heaviside(t-1/2*tao)-Heaviside(-t+1/2*tao) 【例6.4-3】求的Fourier变换，在此是参数，是时间变量。本例演示：fourier的缺省调用格式的使用要十分谨慎；在被变换函数中包含多个符号变量的情况下，对被变换的自变量给予指明，可保证计算结果的正确。syms t x w;ft=exp(-(t-x)*sym(Heaviside(t-x);F1=simple(fourier(ft,t,w)F2=simple(fourier(ft)F3=simple(fourier(ft,t) F1 =1/exp(i*x*w)/(1+i*w)F2 =i*exp(-i*t*w)/(i+w)F3 =i*exp(-t*(2+i*t)/(i+t) 6.4.5.2 Laplace变换及其反变换【例6.4-4】求的Laplace变换。syms t s;syms a b positiveDt=sym(Dirac(t-a);Ut=sym(Heaviside(t-b);Mt=Dt,Ut;exp(-a*t)*sin(b*t),t2*exp(-t);MS=laplace(Mt,t,s) MS = exp(-a*s), exp(-b*s)/s b/(s+a)2+b2), 2/(1+s)3 【例6.4-5】验证Laplace时移性质： 。syms t s;t0=sym(t0,positive);ft=sym(f(t-t0)*sym(Heaviside(t-t0)FS=laplace(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) 6.4.5.3 Z变换及其反变换【例6.4-6】求序列 的Z变换，并用反变换验算。syms nDelta=sym(charfcn0(n);D0=subs(Delta,n,0);D15=subs(Delta,n,15);disp(D0,D15);disp(D0,D15) D0,D15 1, 0 syms z;fn=2*Delta+6*(1-(1/2)n);FZ=simple(ztrans(fn,n,z);disp(FZ = );pretty(FZ),FZ_n=iztrans(FZ,z,n) FZ = 2 4 z + 2 - 2 2 z - 3 z + 1FZ_n =2*charfcn0(n)+6-6*(1/2)n 65 符号代数方程的求解651 线性方程组的符号解【例 6.5-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 X1 = 1 8 8 9 652 一般代数方程组的解S=solve(eq1,eq2,eqn,v1,v2,vn)【例6.5-2】求方程组，关于的解。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 sym z: 2x1 symS.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*w)(1/2)-wS.z 1/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) 【例6.5-3】solve指令求，构成的“欠定”方程组解。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);S.d,S.n,S.p,S.q Warning: 3 equations in 4 variables. In E:MAT53toolboxsymbolicsolve.m at line 110 In E:MAT53toolboxsymbolicsymsolve.m at line 49ans =dans =8ans =4*d+4ans =3*d+6 【例6.5-4】求的解。clear all,syms x;s=solve(x+2)x=2,x) s =.69829942170241042826920133106081 6.6 符号微分方程的求解求微分方程符号解的一般指令S=dsolve(a_1,a_2,a_n)6.6.1 微分方程符号解示例【例6.6-1】求的解。S=dsolve(Dx=y,Dy=-x);disp(blanks(12),x,blanks(21),y),disp(S.x,S.y) x y cos(t)*C1+sin(t)*C2, -sin(t)*C1+cos(t)*C2 【例6.6-2】图示微分方程的通解和奇解的关系。y=dsolve(y=x*Dy-(Dy)2,x)clf,hold on,ezplot(y(2),-6,6,-4,8,1)cc=get(gca,Children);set(cc,Color,r,LineWidth,5)for k=-2:0.5:2;ezplot(subs(y(1),C1,k),-6,6,-4,8,1);endhold off,title(fontname隶书fontsize16通解和奇解) y = x*C1-C12 1/4*x2图 6.6-1 通解和奇解曲线【例6.6-3】求解两点边值问题： 。（注意：相应的数值解法比较复杂）。y=dsolve(x*D2y-3*Dy=x2,y(1)=0,y(5)=0,x) y =-1/3*x3+125/468+31/468*x4 6.7 利用MAPLE的深层符号计算资源6.7.1 MAPLE库函数在线帮助的检索树mhelp index%查看Maple在线帮助的索引类目 Index of help descriptionsCalling Sequence: ?indexcategory or help(index, category);Description:- The following categories of topics are available in the help subsystem:indexexpression#expression operators for forming expressions indexfunction#function list of Maple functions indexmisc#misc miscellaneous facilities indexpackage#packages descriptions of library packages indexprocedure#procedure topics related to procedures and programming indexstatement#statement list of Maple statements To access these help pages, you must prefix the category with index, thus ? indexcategory. mhelp indexfunction Index of descriptions for standard library functionsDescription:- The following are the names of Maples standard library functions. For more information, see ?f where f is any of these