實驗四MATLAB符號計算.doc

實驗四 符號計算符號計算的特點(diǎn)：一，運(yùn)算以推理解析的方式進(jìn)行，因此不受計算誤差積累問題困擾；二，符號計算，或給出完全正確的封閉解，或給出任意精度的數(shù)值解（當(dāng)封閉解不存在時）；三，符號計算指令的調(diào)用比較簡單，經(jīng)典教科書公式相近；四，計算所需時間較長，有時難以忍受。在MATLAB中，符號計算雖以數(shù)值計算的補(bǔ)充身份出現(xiàn)，但涉及符號計算的指令使用、運(yùn)算符操作、計算結(jié)果可視化、程序編制以及在線幫助系統(tǒng)都是十分完整、便捷的。MATLAB的升級和符號計算內(nèi)核Maple的升級，決定著符號計算工具包的升級。但從用戶使用角度看，這些升級所引起的變化相當(dāng)細(xì)微。即使這樣，本章還是及時作了相應(yīng)的更新和說明。如MATLAB 6.5+ 版開始啟用Maple VIII的計算引擎，從而克服了Maple V計算“廣義Fourier變換”時的錯誤（詳見第5.4.1節(jié)）。5.1 符號對象和符號表達(dá)式5.1.1 符號對象的生成和使用【例5.1.1-1】符號常數(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)換為符號變量y=sym(2*sin(x)*cos(x)y=simple(y) y =2*sin(x)*cos(x) y = sin(2*x)【例5.1.1-4】用符號計算驗證三角等式。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】驗證積分。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 符號計算中的算符和基本函數(shù)5.1.3 識別對象類別的指令【例5.1.3-1】數(shù)據(jù)對象及其識別指令的使用。（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 符號表達(dá)式中自由變量的確定【例5.1.4-1】對獨(dú)立自由符號變量的自動辨認(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確定自由變量是對整個矩陣進(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 符號表達(dá)式和符號函數(shù)的操作5.2.1 符號表達(dá)式的操作【例5.2.1-1】按不同的方式合并同冪項。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】對多項式進(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】寫出矩陣各元素的分子、分母多項式（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】簡化（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】簡化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 符號函數(shù)的求反和復(fù)合【例5.2.2-1】求的反函數(shù)syms x;f=x2;g=finverse(f) g =x(1/2) fg=simple(compose(g,f)%驗算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 自動執(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】用簡單算例演示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 符號數(shù)值精度控制和任意精度計算5.2.4.1 向雙精度數(shù)值轉(zhuǎn)換的doblue指令5.2.4.2 任意精度的符號數(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 符號對象與其它數(shù)據(jù)對象間的轉(zhuǎn)換【例5.2.5-1】符號、數(shù)值間的轉(zhuǎn)換。phi=sym(1+sqrt(5)/2)double(phi) phi = 7286977268806824*2(-52) ans = 1.6180【例5.2.5-2】各種多項式表示形式之間的轉(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 符號微積分5.3.1 符號序列的求和【例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 符號微分和矩陣【例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 符號積分5.3.3.1 通用積分指令5.3.3.2 交互式近似積分指令5.3.3.3 符號積分示例【例5.3.3.3-1】求。演示：積分指令對符號函數(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 符號卷積【例5.3.4-1】本例演示卷積的時域積分法：已知系統(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)、輸入同上例，求輸出。對式(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 符號積分變換5.4.1 Fourier變換及其反變換【例 5.4.1-1】求的Fourier變換。本例演示三個重要內(nèi)容：單位階躍函數(shù)和單位脈沖函數(shù)的符號表示；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)行驗算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)%計算結(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ù)，是時間變量。本例演示：fourier的缺省調(diào)用格式的使用要十分謹(jǐn)慎；在被變換函數(shù)中包含多個符號變量的情況下，對被變換的自變量給予指明，可保證計算結(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】驗證Laplace時移性質(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變換，并用反變換驗算。（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 符號代數(shù)方程的求解5.5.1 線性方程組的符號解【例 5.5.1-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 【例 5.5.1-2】求解上例前3 個方程所構(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 符號微分方程的求解5.6.1 符號解法和數(shù)值解法的互補(bǔ)作用5.6.2 求微分方程符號解的一般指令5.6.3 微分方程符號解示例【例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