MATLAB實(shí)用技巧精講之二

第二章MATLAB的數(shù)值計(jì)算功能多項(xiàng)式的運(yùn)算。MATLAB矢量化。MATLAB的符號(hào)計(jì)算功能。1MATLAB的數(shù)值計(jì)算功能：多項(xiàng)式運(yùn)算多項(xiàng)式的表示方法例如:輸入多項(xiàng)式p(x)=x3-5x2+6x-33>>p=[1-56-33];>>poly2sym(p)ans=x^3-5*x^2+6*x-332MATLAB的數(shù)值計(jì)算功能：多項(xiàng)式運(yùn)算求多項(xiàng)式的根用函數(shù)roots求解多項(xiàng)式的根例:求解方程p(x)=2x4-5x3+6x2-x+9的所有根>>p=[2-56-19];>>roots(p)ans=1.6024+1.2709i1.6024-1.2709i-0.3524+0.9755i-0.3524-0.9755i3MATLAB的數(shù)值計(jì)算功能：多項(xiàng)式運(yùn)算多項(xiàng)式的四則運(yùn)算（1）加法和減法

對(duì)于多項(xiàng)式的加法和減法，MATLAB不提供專用的函數(shù)。（2）乘法和除法乘法由函數(shù)conv來實(shí)現(xiàn)，其使用格式為

c=conv(a,b)除法由deconv實(shí)現(xiàn),其使用格式為

d=deconv(a,b)4MATLAB的數(shù)值計(jì)算功能：多項(xiàng)式運(yùn)算多項(xiàng)式的乘除運(yùn)算乘法由函數(shù)conv來實(shí)現(xiàn)除法由deconv實(shí)現(xiàn)例：計(jì)算兩多項(xiàng)式的乘除法>>pd=[6-195432-4539-792-162];>>poly2sym(pd)ans=6*x^6-195*x^5+432*x^4-453*x^3+9*x^2-792*x-162>>d=[3-90-18];>>poly2sym(d)ans=3*x^2-90*x-18>>p1=deconv(pd,d)p1=2-56-195MATLAB的數(shù)值計(jì)算功能：多項(xiàng)式運(yùn)算對(duì)多項(xiàng)式求導(dǎo)數(shù)采用函數(shù)polyder對(duì)多項(xiàng)式求導(dǎo)數(shù)，polyder(p)例：>>p=[6-195432-4539-792-162];>>poly2sym(p)ans=6*x^6-195*x^5+432*x^4-453*x^3+9*x^2-792*x-162>>Dp=polyder(p)Dp=36-9751728-135918-792>>poly2sym(Dp)ans=36*x^5-975*x^4+1728*x^3-1359*x^2+18*x-7926MATLAB的數(shù)值計(jì)算功能：多項(xiàng)式運(yùn)算多項(xiàng)式擬合多項(xiàng)式擬合用擬合函數(shù)polyfit實(shí)現(xiàn)，調(diào)用格式polyfit(X,Y,n)X和Y為擬合數(shù)據(jù)，n為擬合多項(xiàng)式的階數(shù)例：用5階多項(xiàng)式對(duì)[0,π/2]上的正弦函數(shù)進(jìn)行擬合，并圖示擬合情況。>>x=0:pi/20:pi/2;y=sin(x);>>a=polyfit(x,y,5);>>x1=0:pi/30:pi*2;y1=sin(x1);>>y2=a(1)*x1.^5+a(2)*x1.^4+a(3)*x1.^3+a(4)*x1.^2+a(5)*x1+a(6);>>plot(x1,y1,'b-',x1,y2,'r*');>>legend('原曲線','擬合曲線')>>axis([0,7,-1.2,4])7MATLAB的數(shù)值計(jì)算功能：保存與再用把Matlab工作空間中一些有用的數(shù)據(jù)長(zhǎng)久保存下來的方法是生成mat數(shù)據(jù)文件

save——將工作空間中所有的變量存到matlab.mat文件中

savedata將工作空間中所有的變量存到data.mat文件中

savedataab將工作空間中a和b變量存到data.mat文件中

loaddata——將data.mat記錄的變量調(diào)入工作空間8MATLAB的數(shù)值計(jì)算功能：Matlab矢量化何為矢量化問題盡量使用矩陣表述避免出現(xiàn)太多（兩重或以上）的循環(huán)嵌套矢量化技術(shù)的目的在于改善程序的性能例1：有一200×400矩陣datafori=1:200,forj=1:400,ifdata(i,j)>0

data(i,j)=0;endendendcputime=0.0470data(data>0)=0;cputime=0.01609MATLAB的數(shù)值計(jì)算功能：Matlab矢量化例2：上萬個(gè)點(diǎn)的計(jì)算采用矢量化技術(shù)可以大幅度提高運(yùn)算速度%一般循環(huán)編程i=0;fort=0:.01:10000,i=i+1;

y(i)=sin(t);end%矢量化編程t=0:.01:10000;y=sin(t);cputime>20秒cputime=0.2340秒10MATLAB的數(shù)值計(jì)算功能：Matlab矢量化使用向量作索引如果X和V都為向量X(V)就是[X(V(1)),X(V(2)),...,X(V(n))]例：X=[25811141720232629]V=[426]X(V)就是[X(4)X(2)X(6)]，即[11517]11MATLAB的數(shù)值計(jì)算功能：Matlab矢量化創(chuàng)建和操作矩陣從向量構(gòu)建矩陣?yán)?：A=ones(5,5)*10;A=10;A=A(ones(5,5));例2：M=repmat([1:5]',[1,3]);M=

111

222

333

444

555怎樣用repmat生成5×5全10矩陣？乘法是不必要的A=

1010

10

10

10

10

10

10

10

10

10

10

10

10

10

10

10

10

10

10

10

10

10

10

10

12MATLAB的數(shù)值計(jì)算功能：Matlab矢量化排序、設(shè)置和計(jì)數(shù)max最大元素min最小元素sort遞增排序unique遞增排序并且去掉相同元素例1：找出向量中最大元，返回其下標(biāo)>>a=[5278912752578];>>[Y,I]=max(a)Y=525I=713MATLAB的數(shù)值計(jì)算功能：Matlab矢量化排序、設(shè)置和計(jì)數(shù)例2：向量元素排序>>a=[5278912752578];>>sort(a)ans=2577127889525>>unique(a)ans=25712788952514MATLAB的符號(hào)計(jì)算功能1993年之前的符號(hào)計(jì)算語言：Maple、Mathematic和MathCAD等1993年Mathworks公司為了解決MATLAB用戶對(duì)符號(hào)計(jì)算的需求，從加拿大滑鐵盧大學(xué)購(gòu)買了Maple的使用權(quán)，并在此基礎(chǔ)上開發(fā)了符號(hào)計(jì)算工具箱(SymbolicToolbox)至此Matlab

不僅具有數(shù)值運(yùn)算功能，具有符號(hào)計(jì)算功能15MATLAB的符號(hào)計(jì)算功能本部分包含內(nèi)容符號(hào)表達(dá)式、符號(hào)矩陣的創(chuàng)建符號(hào)線性代數(shù)因式分解、展開和簡(jiǎn)化符號(hào)代數(shù)方程求解符號(hào)微積分符號(hào)微分方程16MATLAB的符號(hào)計(jì)算功能

符號(hào)運(yùn)算與數(shù)值運(yùn)算的區(qū)別數(shù)值運(yùn)算中必須先對(duì)變量賦值，然后才能參與運(yùn)算。符號(hào)運(yùn)算無須事先對(duì)獨(dú)立變量賦值，運(yùn)算結(jié)果以標(biāo)準(zhǔn)的符號(hào)形式表達(dá)。符號(hào)運(yùn)算的特點(diǎn)：運(yùn)算對(duì)象可以是沒賦值的符號(hào)變量可以獲得任意精度的解盡管MATLAB7的符號(hào)運(yùn)算功能很強(qiáng)大，但是一門優(yōu)秀的語言應(yīng)該能夠博采眾家之長(zhǎng)。因此MATLAB7語言提供了和MAPLE語言的良好接口，通過maple.m和map.m兩個(gè)專用的M文件來實(shí)現(xiàn)。

17MATLAB符號(hào)變量和符號(hào)表達(dá)式的生成和使用定義基本符號(hào)對(duì)象的指令有兩個(gè)：sym,syms。它們的常用使用格式如下：f=sym(arg)

%把數(shù)字、字符串或表達(dá)式arg

轉(zhuǎn)換為符號(hào)對(duì)象f。syms(‘a(chǎn)rgv1’,‘a(chǎn)rgv2’,‘a(chǎn)rgvk’)%把字符arg1,arg2,argk定義為基本符號(hào)對(duì)象。symsargv1argv2argvk%上述格式的簡(jiǎn)潔形式。18MATLAB符號(hào)變量和符號(hào)表達(dá)式的生成和使用【例】符號(hào)常數(shù)形成中的差異。>>a1=[1/3,pi/7,sqrt(5),pi+sqrt(5)]a1=0.33330.44882.23615.3777>>a2=sym('[1/3,pi/7,sqrt(5),pi+sqrt(5)]')a2=[1/3,pi/7,sqrt(5),pi+sqrt(5)]>>a3=sym([1/3,pi/7,sqrt(5),pi+sqrt(5)])a3=[1/3,pi/7,sqrt(5),6054707603575008*2^(-50)]>>a23=a2-a3a23=[0,0,0,pi+5^(1/2)-189209612611719/35184372088832]19MATLAB符號(hào)變量和符號(hào)表達(dá)式的生成和使用【例】使用sym函數(shù)定義符號(hào)表達(dá)式>>a=sym('a');b=sym('b');c=sym('c');x=sym('x');>>f=a*x^2+b*x+cf=a*x^2+b*x+c也可以采取整體定義法>>f=sym('a*x^2+b*x+c')f=a*x^2+b*x+c20MATLAB符號(hào)變量的生成和使用【例題】用符號(hào)計(jì)算驗(yàn)證三角等式

>>symsfai1fai2>>y=simple(sin(fai1)*cos(fai2)-cos(fai1)*sin(fai2))y=sin(fai1-fai2)>>clear>>y=simple(sin(fai1)*cos(fai2)-cos(fai1)*sin(fai2))???Undefinedfunctionorvariable'fai1'.21符號(hào)表達(dá)式（符號(hào)函數(shù)）的操作符號(hào)表達(dá)式的四則運(yùn)算>>symsxyab>>fun1=sin(x)+cos(y);>>fun2=a+b;>>fun1+fun2ans=sin(x)+cos(y)+a+b>>fun1*fun2ans=(sin(x)+cos(y))*(a+b)>>fun1/fun2ans=(sin(x)+cos(y))/(a+b)22

MATLAB符號(hào)矩陣的創(chuàng)建

用Matlab函數(shù)sym創(chuàng)建矩陣命令格式：A=sym('[]')

※需用sym函數(shù)※需用''標(biāo)識(shí)注意：數(shù)值矩陣A=[1,2;3,4]有效，A=[a,b;c,d]

則不識(shí)別！例：>>B=sym('[a,2*b;3*a,0]')B=

[a,2*b]

[3*a,0]這就完成了一個(gè)符號(hào)矩陣的創(chuàng)建。符號(hào)矩陣的每一行的兩端都有方括號(hào)，這是與Matlab數(shù)值矩陣的一個(gè)重要區(qū)別。23MATLAB符號(hào)矩陣的創(chuàng)建用生成子矩陣的方法生成符號(hào)矩陣命令格式：A=[‘[]’；‘[]’]

【例題】用生成子矩陣的方法生成符號(hào)矩陣>>h1=['[55,ttt]';'[44,1]']???Errorusing==>vertcatAllrowsinthebracketedexpressionmusthavethesamenumberofcolumns.>>h2=['[55,ttt]';'[44,1]']h2=[55,ttt][44,1]該方法不需要調(diào)用sym函數(shù)，但是要保證同一列的元素具有相同的長(zhǎng)度。當(dāng)要輸入的字符串長(zhǎng)度不一致時(shí)，可以在較短的字符串前邊或者后邊加上空格符補(bǔ)充24MATLAB符號(hào)矩陣的創(chuàng)建由數(shù)值矩陣轉(zhuǎn)換成符號(hào)矩陣將數(shù)值矩陣轉(zhuǎn)化為符號(hào)矩陣，函數(shù)調(diào)用格式：B=sym(A)例：>>A=[1/3,2.5;1/0.7,2/5]A=0.33332.50001.42860.4000>>B=sym(A)B=[1/3,5/2][10/7,2/5]>>A+B=?ans=[2/3,5][20/7,4/5]25MATLAB符號(hào)矩陣及符號(hào)數(shù)組的運(yùn)算符號(hào)矩陣的四則運(yùn)算A+B和A-B命令可以實(shí)現(xiàn)符號(hào)陣列的加法和減法。A*B命令可以實(shí)現(xiàn)符號(hào)矩陣的乘法。A\B命令可以實(shí)現(xiàn)矩陣的左除法。

A\B為符號(hào)線性方程組A*X=B的解。B/A命令可以實(shí)現(xiàn)矩陣的右除法。

B/A為符號(hào)線性方程組X*A=B的解。26MATLAB符號(hào)矩陣及符號(hào)數(shù)組的運(yùn)算【例】>>m=sym('[a,a^2,a*2,1/a]')m=[a,a^2,a*2,1/a]>>n=sym('[2*a,b,a,a^2]')n=[2*a,b,a,a^2]>>m+nans=[3*a,a^2+b,3*a,1/a+a^2]27MATLAB符號(hào)矩陣及符號(hào)數(shù)組的運(yùn)算>>x=m\nx=[0,0,0,0][0,0,0,0][0,0,0,0][2*a^2,b*a,a^2,a^3]>>y=m/nInfWarning:Systemisinconsistent.Solutiondoesnotexist.y=Inf28MATLAB符號(hào)矩陣及符號(hào)數(shù)組的運(yùn)算符號(hào)矩陣的逆和行列式運(yùn)算

inv函數(shù)可以用來求方正的逆矩陣，命令為：inv(X)；

det函數(shù)可以求方正的行列式，命令為：det(X)。

29MATLAB符號(hào)矩陣及符號(hào)數(shù)組的運(yùn)算符號(hào)數(shù)組的四則運(yùn)算A.*B命令用于符號(hào)數(shù)組的乘法。A.*B為按參量A與B對(duì)應(yīng)的分量進(jìn)行相乘。A.\B命令用于數(shù)組的左除法運(yùn)算。A.\B為按對(duì)應(yīng)的分量進(jìn)行相除。A./B命令用于數(shù)組的右除法運(yùn)算。A./B為按對(duì)應(yīng)的分量進(jìn)行相除。30MATLAB符號(hào)矩陣及符號(hào)數(shù)組的運(yùn)算【例】符號(hào)數(shù)組的四則運(yùn)算>>q=sym('[3,4,9,6;x,y,z,w;a,b,c,d]')

q=[3,4,9,6][x,y,z,w][a,b,c,d]>>p=sym('[x,1/x,x^2,x^3;a,b,c,d;5,2,3,6]')p=[x,1/x,x^2,x^3][a,b,c,d][5,2,3,6]>>r=q*p???Errorusing==>sym.mtimesInnermatrixdimensionsmustagree.31MATLAB符號(hào)矩陣及符號(hào)數(shù)組的運(yùn)算>>q.*pans=[3*x,4/x,9*x^2,6*x^3][x*a,y*b,z*c,w*d][5*a,2*b,3*c,6*d]>>q.\p

ans=[1/3*x,1/4/x,1/9*x^2,1/6*x^3][a/x,b/y,c/z,d/w][5/a,2/b,3/c,6/d]>>q./pans=[3/x,4*x,9/x^2,6/x^3][x/a,y/b,z/c,w/d][1/5*a,1/2*b,1/3*c,1/6*d]32MATLAB符號(hào)極限符號(hào)極限limit(F,a)命令中設(shè)定變量為x，計(jì)算當(dāng)時(shí)F的極限值。limit(F)命令中設(shè)定變量為x，計(jì)算當(dāng)時(shí)F的極限值。limit(F,x,a,`right`)或limit(F,x,a,`left`)來計(jì)算符號(hào)函數(shù)F的單側(cè)極限：左極限或右極限。33MATLAB符號(hào)極限【例】使用limit函數(shù)來求符號(hào)極限。>>symsxath;>>limit(sin(x)/x)ans=1>>limit((x-2)/(x^2-4),2)ans=1/4>>limit(1/x,x,0,'right')ans=Inf>>limit(1/x,x,0,'left')ans=-Inf34MATLAB符號(hào)微積分符號(hào)微積分diff(f)——對(duì)缺省變量求一階導(dǎo)數(shù)diff(f,v)——對(duì)指定變量v求一階導(dǎo)數(shù)diff(f,v,n)——對(duì)指定變量v求n階導(dǎo)數(shù)int(f)——對(duì)f表達(dá)式的缺省變量求積分int(f,v)——對(duì)f表達(dá)式的v變量求積分int(f,v,a,b)——對(duì)f表達(dá)式的v變量在(a,b)區(qū)間求定積分int(‘被積表達(dá)式’，‘積分變量’，‘積分上限’，‘積分下限’)——定積分——缺省時(shí)為不定積分35符號(hào)微分和求導(dǎo)【例】求單個(gè)自變量時(shí)的微分>>symsx>>diff(x^3+2*x^2+4*x+6)

ans=3*x^2+4*x+4

>>diff(sin(x^3),4)

ans=81*sin(x^3)*x^8-324*cos(x^3)*x^5-180*sin(x^3)*x^236符號(hào)微分和求導(dǎo)【例】對(duì)多個(gè)自變量的函數(shù)中的某個(gè)自變量求導(dǎo)。>>symsxy>>diff(x*y+y^2+sin(x)+cos(y),y)ans=x+2*y-sin(y)>>diff(x*y+y^2+sin(x)+cos(y),y,3)ans=sin(y)37MATLAB的符號(hào)積分例：計(jì)算二重不定積分38MATLAB的符號(hào)積分例：計(jì)算二重不定積分>>symsxy>>F=int(int(x*exp(-x*y),x),y)

F=

1/y*exp(-x*y)39MATLAB符號(hào)積分變換Fourier變換及其逆變換

時(shí)域中的f(t)與它在頻域中的Fourier變換F(w)之間存在如下關(guān)系：

Fw=fourier(ft,t,w)%求時(shí)域函數(shù)ft的Fourier變換Fwft=ifourier(Fw,w,t)%求頻域函數(shù)Fw的Fourier反變換ft40MATLAB符號(hào)積分變換Fourier變換及其逆變換【例】求的Fourier變換。>>symstw>>ut=sym('Heaviside(t)')>>UT=fourier(ut)UT=pi*dirac(w)-i/w

>>vt=ifourier(UT,w,t)

vt=heaviside(t)41MATLAB符號(hào)積分變換Laplace變換及其逆變換的定義為：

Fs=laplace(ft,t,s)%求時(shí)域函數(shù)ft的Laplace變換Fsft=ilaplace(Fs,s,t)%求頻域函數(shù)Fs的Laplace反變換ft42MATLAB符號(hào)積分變換【例】求的Laplace變換。>>

symsts;>>

symsabpositive;>>

Dt=sym('Dirac(t-a)');>>

Ut=sym('Heaviside(t-b)');>>Mt=[Dt,Ut;exp(-a*t)*sin(b*t),t^2*exp(-t)];>>Ms=laplace(Mt,t,s)Ms=[exp(-s*a),exp(-s*b)/s][1/b/((s+a)^2/b^2+1),2/(s+1)^3]43

MATLAB符號(hào)代數(shù)方程的求解

符號(hào)代數(shù)方程求解Matlab符號(hào)運(yùn)算能夠解一般的線性方程、非線性方程及一般的代數(shù)方程、代數(shù)方程組。命令格式：solve(f)——求一個(gè)方程的解solve(f1,f2,…fn)——求n個(gè)方程的解solve(‘eq1’,‘eq2’,…,‘eqn’,‘v1’,‘v2’,…,‘vn’)——求方程組關(guān)于指定變量的解（推薦格式）。solve(exp1,exp2,….,expn,v1,v2,…,vn)——求方程組關(guān)于指定變量的解（推薦格式)。44

MATLAB符號(hào)代數(shù)方程的求解

符號(hào)代數(shù)方程求解例1：f=ax2+bx+c求解>>f='a*x^2+b*x+c';>>solve(f)%對(duì)缺省變量x求解ans=[1/2/a*(-b+(b^2-4*a*c)^(1/2))][1/2/a*(-b-(b^2-4*a*c)^(1/2))]計(jì)算機(jī)格式一般格式45MATLAB符號(hào)代數(shù)方程的求解符號(hào)代數(shù)方程求解>>solve(f,'b')%對(duì)指定變量b求解ans=-(a*x^2+c)/x例2：符號(hào)方程cos(x)=sin(x)求解>>f1=solve('cos(x)=sin(x)')f1=1/4*pi46MATLAB符號(hào)代數(shù)方程的求解符號(hào)代數(shù)方程求解例：解方程組>>g1='x+y+z=1';g2='x-y+z=2';g3='2*x-y-z=1';>>f=solve(g1,g2,g3)f=x:[1x1sym]

y:[1x1sym]z:[1x1sym]>>[f.x,f.y,f.z]ans=[2/3,-1/2,5/6]f=solve('x+y+z=1','x-y+z=2','2*x-y-z=1')[x,y,z]=solve('x+y+z=1','x-y+z=2','2*x-y-z=1')x=2/3y=-1/2z=5/647MATLAB符號(hào)微分方程求解符號(hào)微分方程求解用一個(gè)函數(shù)可以方便地得到微分方程的符號(hào)解，求解指令：dsolve調(diào)用格式：dsolve(‘f’,‘g’)f微分方程，可多至12個(gè)微分方程的求解；g為初始條件默認(rèn)自變量為'x',可任意指定