實(shí)驗(yàn)四MATLAB符號(hào)計(jì)算.doc

實(shí)驗(yàn)四 符號(hào)計(jì)算符號(hào)計(jì)算的特點(diǎn)：一，運(yùn)算以推理解析的方式進(jìn)行，因此不受計(jì)算誤差積累問題困擾；二，符號(hào)計(jì)算，或給出完全正確的封閉解，或給出任意精度的數(shù)值解（當(dāng)封閉解不存在時(shí)）；三，符號(hào)計(jì)算指令的調(diào)用比較簡(jiǎn)單，經(jīng)典教科書公式相近；四，計(jì)算所需時(shí)間較長(zhǎng)，有時(shí)難以忍受。在MATLAB中，符號(hào)計(jì)算雖以數(shù)值計(jì)算的補(bǔ)充身份出現(xiàn)，但涉及符號(hào)計(jì)算的指令使用、運(yùn)算符操作、計(jì)算結(jié)果可視化、程序編制以及在線幫助系統(tǒng)都是十分完整、便捷的。MATLAB的升級(jí)和符號(hào)計(jì)算內(nèi)核Maple的升級(jí)，決定著符號(hào)計(jì)算工具包的升級(jí)。但從用戶使用角度看，這些升級(jí)所引起的變化相當(dāng)細(xì)微。即使這樣，本章還是及時(shí)作了相應(yīng)的更新和說明。如MATLAB 6.5+ 版開始啟用Maple VIII的計(jì)算引擎，從而克服了Maple V計(jì)算“廣義Fourier變換”時(shí)的錯(cuò)誤（詳見第5.4.1節(jié)）。5.1 符號(hào)對(duì)象和符號(hào)表達(dá)式5.1.1 符號(hào)對(duì)象的生成和使用【例5.1.1-1】符號(hào)常數(shù)形成中的差異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),e)%a4=sym(1/3,pi/7,sqrt(5),pi+sqrt(5)%a24=a2-a4 a1 = 0.3333 0.4488 2.2361 5.3777a2 = 1/3, pi/7, sqrt(5), 6054707603575008*2(-50)a3 = 1/3-eps/12, pi/7-13*eps/165, sqrt(5)+137*eps/280, 6054707603575008*2(-50) a4 = 1/3, pi/7, sqrt(5), pi+sqrt(5) a24 = 0, 0, 0, 189209612611719/35184372088832-pi-5(1/2) 【例5.1.1-2】演示：幾種輸入下產(chǎn)生矩陣的異同。a1=sym(1/3,0.2+sqrt(2),pi)%a2=sym(1/3,0.2+sqrt(2),pi)%a3=sym(1/3 0.2+sqrt(2) pi)%a1_a2=a1-a2% a1 = 1/3, 7269771597999872*2(-52), pia2 = 1/3, 0.2+sqrt(2), pi a3 = 1/3, 0.2+sqrt(2), pi a1_a2 = 0, 1.4142135623730951010657008737326-2(1/2), 0 【例5.1.1-3】把字符表達(dá)式轉(zhuǎn)換為符號(hào)變量y=sym(2*sin(x)*cos(x)y=simple(y) y =2*sin(x)*cos(x) y = sin(2*x)【例5.1.1-4】用符號(hào)計(jì)算驗(yàn)證三角等式。syms fai1 fai2;y=simple(sin(fai1)*cos(fai2)-cos(fai1)*sin(fai2) y =sin(fai1-fai2)【例5.1.1-5】求矩陣的行列式值、逆和特征根syms a11 a12 a21 a22;A=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) 【例5.1.1-6】驗(yàn)證積分。syms A t tao w;yf=int(A*exp(-i*w*t),t,-tao/2,tao/2);Yf=simple(yf) Yf =2*A*sin(1/2*w*tao)/w 5.1.2 符號(hào)計(jì)算中的算符和基本函數(shù)5.1.3 識(shí)別對(duì)象類別的指令【例5.1.3-1】數(shù)據(jù)對(duì)象及其識(shí)別指令的使用。（1）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（2）SizeMn=size(Mn),SizeMc=size(Mc),SizeMs=size(Ms) SizeMn = 2 2SizeMc = 1 9SizeMs = 2 2（3）CMn=class(Mn),CMc=class(Mc),CMs=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 Mc 1x9 18 char array Mn 2x2 32 double array Ms 2x2 312 sym objectGrand total is 21 elements using 362 bytes5.1.4 符號(hào)表達(dá)式中自由變量的確定【例5.1.4-1】對(duì)獨(dú)立自由符號(hào)變量的自動(dòng)辨認(rèn)。（1）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; （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【例5.1.4-2】findsym確定自由變量是對(duì)整個(gè)矩陣進(jìn)行的。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)+v ans =x5.2 符號(hào)表達(dá)式和符號(hào)函數(shù)的操作5.2.1 符號(hào)表達(dá)式的操作【例5.2.1-1】按不同的方式合并同冪項(xiàng)。EXPR=sym(x2+x*exp(-t)+1)*(x+exp(-t);expr1=collect(EXPR)expr2=collect(EXPR,exp(-t) expr1 = x3+2*exp(-t)*x2+(1+exp(-t)2)*x+exp(-t) expr2 = x*exp(-t)2+(2*x2+1)*exp(-t)+(x2+1)*x【例5.2.1-2】factor指令的使用（1）syms a x;f1=x4-5*x3+5*x2+5*x-6;factor(f1) ans = (x-1)*(x-2)*(x-3)*(x+1)（2）f2=x2-a2;factor(f2) ans = -(a-x)*(a+x)（3）factor(1025) ans = 5 5 41【例5.2.1-3】對(duì)多項(xiàng)式進(jìn)行嵌套型分解clear;syms a x;f1=x4-5*x3+5*x2+5*x-6;horner(f1) ans = -6+(5+(5+(-5+x)*x)*x)*x【例5.2.1-4】寫出矩陣各元素的分子、分母多項(xiàng)式（1）syms x;A=3/2,(x2+3)/(2*x-1)+3*x/(x-1);4/x2,3*x+4;n,d=numden(A)pretty(simplify(A)% n = 3, x3+5*x2-3 4, 3*x+4d = 2, (2*x-1)*(x-1) x2, 1 3 2 x + 5 x - 3 3/2 - (2 x - 1) (x - 1) 4 - 3 x + 4 2 x （2）pretty(simplify(n./d) 3 2 x + 5 x - 3 3/2 - (2 x - 1) (x - 1) 4 - 3 x + 4 2 x 【例5.2.1-5】簡(jiǎn)化（1）syms x;f=(1/x3+6/x2+12/x+8)(1/3);sfy1=simplify(f),sfy2=simplify(sfy1) sfy1 = (2*x+1)3/x3)(1/3)sfy2 = (2*x+1)3/x3)(1/3)（2）g1=simple(f),g2=simple(g1) g1 = (2*x+1)/x g2 =2+1/x【例5.2.1-6】簡(jiǎn)化syms x;ff=cos(x)+sqrt(-sin(x)2);ssfy1=simplify(ff),ssfy2=simplify(ssfy1) ssfy1 = cos(x)+(-1+cos(x)2)(1/2) ssfy2 = cos(x)+(-1+cos(x)2)(1/2)gg1=simple(ff),gg2=simple(gg1) gg1 = cos(x)+i*sin(x) gg2 = exp(i*x)5.2.2 符號(hào)函數(shù)的求反和復(fù)合【例5.2.2-1】求的反函數(shù)syms x;f=x2;g=finverse(f) g =x(1/2) fg=simple(compose(g,f)%驗(yàn)算g(f(x)是否等于x fg =x【例5.2.2-2】求的復(fù)合函數(shù)（1）syms x y u fai t;f=x/(1+u2);g=cos(y+fai);fg1=compose(f,g) fg1 =cos(y+fai)/(1+u2)（2）fg2=compose(f,g,u,fai,t) fg2 = x/(1+cos(y+t)2)5.2.3 置換及其應(yīng)用5.2.3.1 自動(dòng)執(zhí)行的子表達(dá)式置換指令【例5.2.3.1-1】演示子表達(dá)式的置換表示。clear all,syms a b c d W;V,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)5.2.3.2 通用置換指令【例5.2.3.2-1】用簡(jiǎn)單算例演示subs的置換規(guī)則。（1）syms a x;f=a*sin(x)+5; f =a*sin(x)+5（2）f1=subs(f,sin(x),sym(y)% f1 = a*y+5（3）f2=subs(f,a,x,2,sym(pi/3)% f2 = 3(1/2)+5（4）f3=subs(f,a,x,2,pi/3)% f3 = 6.7321（5）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（6）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.00005.2.4 符號(hào)數(shù)值精度控制和任意精度計(jì)算5.2.4.1 向雙精度數(shù)值轉(zhuǎn)換的doblue指令5.2.4.2 任意精度的符號(hào)數(shù)值【例5.2.4.2-1】指令使用演示。digits Digits = 32p0=sym(1+sqrt(5)/2); p1=sym(1+sqrt(5)/2)e01=vpa(abs(p0-p1) p1 =7286977268806824*2(-52) e01 =.543211520368251e-16p2=vpa(p0)e02=vpa(abs(p0-p2),64) p2 = 1.6180339887498948482045868343656 e02 = .38117720309179805762862135448622e-31digits Digits = 325.2.5 符號(hào)對(duì)象與其它數(shù)據(jù)對(duì)象間的轉(zhuǎn)換【例5.2.5-1】符號(hào)、數(shù)值間的轉(zhuǎn)換。phi=sym(1+sqrt(5)/2)double(phi) phi = 7286977268806824*2(-52) ans = 1.6180【例5.2.5-2】各種多項(xiàng)式表示形式之間的轉(zhuǎn)換syms x;f=x3+2*x2-3*x+5;sy2p=sym2poly(f)p2st=poly2str(sy2p,x)p2sy=poly2sym(sy2p)pretty(f,x) sy2p = 1 2 -3 5p2st = x3 + 2 x2 - 3 x + 5p2sy = x3+2*x2-3*x+55.3 符號(hào)微積分5.3.1 符號(hào)序列的求和【例5.3.1-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*t s2 = 1/8*pi2, -log(2)5.3.2 符號(hào)微分和矩陣【例5.3.2-1】求、和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/x dfdt2 = 0, 6*t 0, 0 dfdxdt = 0, 0 -sin(x), 0【例5.3.2-2】求的矩陣。syms x1 x2 x3;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)5.3.3 符號(hào)積分5.3.3.1 通用積分指令5.3.3.2 交互式近似積分指令5.3.3.3 符號(hào)積分示例【例5.3.3.3-1】求。演示：積分指令對(duì)符號(hào)函數(shù)矩陣的作用。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) 【例5.3.3.3-2】求。演示如何使用mfun指令獲取一組積分值。（1）F1=int(1/log(t),t,0,x) F1 =-Ei(1,-log(x)（2）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【例5.3.3.3-3】求積分。注意：內(nèi)積分上下限都是函數(shù)。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 = 64/225*2(3/4)-6072064/348075*2(1/2)+14912/4641*2(1/4)+1610027357/6563700 VF2 = 224.92153573331143159790710032805【例5.3.3.3-4】利用rsums求積分。（與例5.3.3.3-2結(jié)果比較）syms x positive;px=0.5/log(0.5*x);rsums(px) 圖5.3-1 5.3.4 符號(hào)卷積【例5.3.4-1】本例演示卷積的時(shí)域積分法：已知系統(tǒng)沖激響應(yīng)，求輸入下的輸出響應(yīng)。syms T t tao;ut=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 = -1/(T-1)/exp(t)+1/(T-1)/exp(t/T)【例5.3.4-2】本例演示通過變換和反變換求取卷積。系統(tǒng)沖激響應(yīng)、輸入同上例，求輸出。對(duì)式(5.3.4-1)兩邊進(jìn)行Laplace變換得，因此有syms s;yt=ilaplace(laplace(ut,t,s)*laplace(ht,t,s),s,t);yt=simple(yt) yt = (exp(-t/T)-exp(-t)/(T-1)【例5.3.4-3】求函數(shù)和的卷積。syms tao;t=sym(t,positive);ut=sym(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 = (exp(1-t)*t-1)*heaviside(t-1)+1+(-t-1)/exp(t)5.4 符號(hào)積分變換5.4.1 Fourier變換及其反變換【例 5.4.1-1】求的Fourier變換。本例演示三個(gè)重要內(nèi)容：?jiǎn)挝浑A躍函數(shù)和單位脈沖函數(shù)的符號(hào)表示；fourier指令的使用；simple指令的表現(xiàn)。（1）求Fourier變換syms t w;ut=sym(Heaviside(t);%UT=fourier(ut)UTC=maple(convert,UT,piecewise,w)%UTS=simple(UT) UT = 2*pi*dirac(w)UTC =PIECEWISE(pi*NaN, w = 0,0, otherwise)UTS = 2*pi*dirac(w)（2）求Fourier反變換進(jìn)行驗(yàn)算Ut=ifourier(UT,w,t)Uts=ifourier(UTS,w,t) Ut = 1Uts =1【例5.4.1-2】用fourier指令求例5.1.1-6中方波脈沖的Fourier變換。本例演示：fourier, simple 指令的配合使用。（1）syms A t wsyms tao positive%yt=sym(Heaviside(t+tao/2)-Heaviside(t-tao/2);Yw=fourier(A*yt,t,w)Ywc=maple(convert,Yw,piecewise,w)%計(jì)算結(jié)果起指示作用Yws=simple(Yw) Yw =A*(i*exp(-1/2*i*tao*w)/w+pi*dirac(w)Ywc = PIECEWISE(A*(i*exp(-1/2*i*tao*w)/w+pi*NaN), w = 0,i*A*exp(-1/2*i*tao*w)/w, otherwise) Yws =A*(i*exp(-1/2*i*tao*w)/w+pi*dirac(w)（2）Yt=ifourier(Yw,w,t)Yst=ifourier(Yws,w,t) Yt = A*heaviside(-t+1/2*tao)Yst = A*heaviside(-t+1/2*tao)【例5.4.1-3】求的Fourier變換，在此是參數(shù)，是時(shí)間變量。本例演示：fourier的缺省調(diào)用格式的使用要十分謹(jǐn)慎；在被變換函數(shù)中包含多個(gè)符號(hào)變量的情況下，對(duì)被變換的自變量給予指明，可保證計(jì)算結(jié)果的正確。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 = fourier(exp(-t),t,w)/exp(i*x*w) F2 =exp(-t)*fourier(exp(x)*heaviside(t-x),x,w) F3 =exp(-t)*fourier(exp(x)*heaviside(t-x),x,t)5.4.2 Laplace變換及其反變換【例5.4.2-1】求的Laplace變換。syms t s;syms a b positive%Dt=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(-s*a), exp(-s*b)/s 1/b/(s+a)2/b2+1), 2/(1+s)3【例5.4.2-2】驗(yàn)證Laplace時(shí)移性質(zhì)： 。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)5.4.3 Z變換及其反變換【例5.4.3-1】求序列 的Z變換，并用反變換驗(yàn)算。（1）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（2）求序列的Z變換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 + 1 FZ_n =2*charfcn0(n)+6-6*(1/2)n5.5 符號(hào)代數(shù)方程的求解5.5.1 線性方程組的符號(hào)解【例 5.5.1-1】求線性方程組的解。本例演示，符號(hào)線性方程組的基本解法。該方程組的矩陣形式是 。該式簡(jiǎn)記為。求符號(hào)解的指令如下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 【例 5.5.1-2】求解上例前3 個(gè)方程所構(gòu)成的“欠定”方程組，并解釋解的含義。A2=A(1:3,:);X2=A2b(1:3,1)syms k;XX2=X2+k*null(A2)A2*XX2 X1 =18895.5.2 一般代數(shù)方程組的解【例5.5.2-1】求方程組，關(guān)于的解。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)-w S.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)【例5.5.2-2】用solve指令重做例 5.5.1-2。即求，構(gòu)成的“欠定”方程組解。syms d n p q;eq1=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 ans = d ans = 8ans =4*d+4 ans = 3*d+6【例5.5.2-3】求的解。clear all,syms x;s=solve(x+2)x=2,x) s =.698299421702410428269201331060815.6 符號(hào)微分方程的求解5.6.1 符號(hào)解法和數(shù)值解法的互補(bǔ)作用5.6.2 求微分方程符號(hào)解的一般指令5.6.3 微分方程符號(hào)解示例【例5.6.3-1】求的解。S=dsolve(Dx=y,Dy=-x);disp(blanks(12),x,blanks(21),y),disp(S.x,S.y) x y C1*sin(t)+C2*cos(t), C1*cos(t)-C2*sin(t)【例5.6.3-2】圖示微分方程的通解和奇解的關(guān)系。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通解和奇解) 圖 5.6-1 【例5.6.3-3】求解兩點(diǎn)邊值問題： 。（注意：相應(yīng)的數(shù)值解法比較復(fù)雜）。y=dsolve(x*D2y-3*Dy=x2,y(1)=0,y(5)=0,x) y = 31/468*x4-1/3*x3+125/468 【例5.6.3-4】求邊值問題的解。（注意：相應(yīng)的數(shù)值解法比較復(fù)雜）。S=dsolve(Df=3*f+4*g,Dg=-4*f+3*g,f(0)=0,f(3)=1)S.f,S.g S = f: 1x1 sym g: 1x1 symans =exp(3*t)/sin(12)/(cosh(9)+si