第5章MATLAB程序設計教學教材

第5章

MATLAB程序設計

MATLAB有兩種常用工作方式：（1)直接交互的指令行操作方式；（2)M文件中的編程工作方式。5.1M文件中的功能和特點

(1)MATLAB程序文件中一個ASCⅡ碼文件，擴展名一律為“.m”，可用任何字處理軟件進行編輯和修改。

(2)

MATLAB 是程序性的語言是解釋性語言。

(3)M文件大大擴展了MATLAB的能力5.2M文件的建立與使用1、M文件的形式

M文件中有兩種形式

(1)命令文件

（ScriptFile）

也叫做腳本文件

(2)函數(shù)文件（FunctionFile）2、腳本M文件

（1）

建立

用記事本（Notebook）編寫程序；

在MATLAB系統(tǒng)下選擇：“file”→“new”→“m-file”

（2）內(nèi)容：按照程序的功能，依據(jù)MATLAB的程序結構，組織的合法的MATLAB命令。(4)注意

（1）命令文件中的語句可以訪問MATLAB工作區(qū)(workspace)中的所有變量。運行過程中所有產(chǎn)生的變量均是全局變量。（2）“%”開始的行為注釋行，不予執(zhí)行。（3）若將文件未存放在系統(tǒng)的搜索目錄下時，運用該文件之前，應先進入該目錄。（如cdd:\mywok）(3)運行

在命令窗口直接鍵入文件名.m3、函數(shù)文件(1)建立 用記事本（Notebook）編寫程序 在MATLAB系統(tǒng)下選擇：“file”→“new”→“m-file”(2)函數(shù)文件的格式Function輸出參數(shù)=函數(shù)名（輸入?yún)?shù)表）函數(shù)體(3)函數(shù)體的內(nèi)容

函數(shù)定義行

Hi行：幫助信息的第一行，用于提示函數(shù)的功能，當用命令Lookfor查詢該函數(shù)的幫助信息時將顯示該行內(nèi)容。

幫助體

函數(shù)體函數(shù)名命名規(guī)則同變量名。以字母、下劃線和數(shù)字組成，不識別函數(shù)名共31個字符。輸入?yún)?shù)：用中括號“[]”括起來，參數(shù)兩兩之間用逗號隔開。輸出參數(shù)：函數(shù)m文件的運算結果傳到調(diào)用處，當參數(shù)不只一個時，用逗號隔開。函數(shù)文件名由函數(shù)名再加后綴“m”組成。函數(shù)文件中定義變量為局部變量，函數(shù)文件命令運行結束，該類變量會自動釋放。如果變量在程序運行前就已存在的話，程序運行后它不會受到影響(4)注意

（5）應用舉例函數(shù)表示函數(shù)m文件的內(nèi)容為：function

y=foft(x)y=1./(1+exp(-x))文件名為：foftfplot(＇foft＇,[020],1e-04)數(shù)學函數(shù)曲線的繪圖fplot(fun，lims)其中：fun表示函數(shù)名，定義函數(shù)的M文件名lims=[XMINXMAXYMINYMAX]函數(shù)的極小值點和零值點l

單變量函數(shù)的局部極小值點

fmin(‘F’，x1，x2，options)F函數(shù)名，x1,x2指定極小值區(qū)間，options第一個分量為正時，則顯示函數(shù)的運行步驟。Options（1）屏幕上是否顯示最小值的迭代過程，非0時，顯示，為0時不顯示，默認為0Options（2）迭代時自變量的誤差控程限，默認為1.0e-04Options（3）為迭代時函數(shù)值的誤差控程限，函數(shù)fmins使用該參數(shù)Options（4）為迭代時最大的迭代次數(shù)，函數(shù)fmin默認值為500，fmins為200步。舉例：定義函數(shù)文件humps.mfunctiony=humps(x)y=1./((x-0.3).^2+0.01)+1./((x-0.9).^2+0.04)-6figure(1)subplot(2,2,1),fplot(‘humps’,[-5,5]),gridonsubplot(2,2,2),fplot(‘humps’,[-5,5-1025]),gridonsubplot(2,2,3),fplot(‘5*sin(x)’,[-1,1]),gridonsubplot(2,2,3),fplot(‘[5*sin(x),humps]’,[-1,1])gridonx=fmin(‘humps’,0.3,1)極小值點x函數(shù)的極小值為：humps（x）

l

函數(shù)的零點fzero(‘F’，猜測值)得到零點的區(qū)間l

多變量函數(shù)的局部極小值點fmins(‘函數(shù)名’，初始極值向量，options)得到向量附近的一個極小值舉例：定義函數(shù)

function

b=two_var(v)x=v(1);y=v(2)b=x.^2+y.^2-x.^2*y.^2給定初始值：v=[-0.7–1.2]調(diào)用函數(shù)：a=fmins(‘two_var’,v)options與單變量時相同

數(shù)值積分

quad(‘F’,A,B,TOL,TRACE)

其中：‘F’為函數(shù)名，A和B指定積分區(qū)間，TOL迭代誤差，TRACE，采用的遞歸的辛浦生方法。quad8(‘F’,A,B,TOL,TRACE)

其中：‘F’為函數(shù)名，A和B指定積分區(qū)間，TOL迭代誤差，TRACE，采用的newton-cotes方法。dblquad(‘F’,inmin,inmax,outnmin,outinmax)其中：

inmin與inmax，分別為內(nèi)層積分的上下限，outnmin與outinmax，分別為外層積分的上下限舉例：定義函數(shù)

functiony=ex542f(x)y=-x.*x+115

調(diào)用函數(shù)：

s=quad(‘ex542f’,0,10)

例建立一個命令文件將變量a,b的值互換，然后運行該命令文件。首先建立命令文件并以文件名exch.m存盤：clear;a=1:10;b=[11,12,13,14;15,16,17,18];c=a;a=b;b=c;ab然后在MATLAB的命令窗口中輸入exch，將會執(zhí)行該命令文件。例建立一個函數(shù)文件將變量a,b的值互換，然后在命令窗口調(diào)用該函數(shù)文件。首先建立函數(shù)文件fexch.m：function[a,b]=exch(a,b)c=a;a=b;b=c;然后在MATLAB的命令窗口調(diào)用該函數(shù)文件：clear;x=1:10;y=[11,12,13,14;15,16,17,18];[x,y]=fexch(x,y)5.2

數(shù)據(jù)的輸入輸出1、input函數(shù)格式：A=input(提示信息，選項)；其中：提示信息為一個字符串，用于提示用戶輸入什么樣的數(shù)據(jù)。如果在input函數(shù)調(diào)用時采用's'選項，則允許用戶輸入一個字符串。例如，想輸入一個人的姓名，可采用命令：xm=input('What''syourname?','s')

2、disp函數(shù)格式:disp(輸出項)其中輸出項既可以為字符串，也可以為矩陣。注意：用disp函數(shù)顯示矩陣時將不顯示矩陣的名字，而且其格式更緊密，且不留任何沒有意義的空行。

例1求一元二次方程ax2+bx+c=0的根。程序如下：a=input('a=?');b=input('b=?');c=input('c=?');d=b*b-4*a*c;x=[(-b+sqrt(d))/(2*a),(-b-sqrt(d))/(2*a)];disp(['x1=',num2str(x(1)),',x2=',num2str(x(2))]);5.3M文件的程序結構1、順序結構2、分支結構（1）if（條件）（命令組）end（2)if（條件）（命令組1）else（命令2）end當條件成立時，則執(zhí)行語句組，執(zhí)行完之后繼續(xù)執(zhí)行if語句的后繼語句，若條件不成立，則直接執(zhí)行if語句的后繼語句

2、分支結構(3)if

elseif結構

if（條件1）命令組1elseif(條件2)命令組2else命令序列end例2

刪除數(shù)據(jù)中的無關數(shù)據(jù)：一般認為，誤差大于３倍標準差的數(shù)據(jù)項，就認為是無關項，因此可以從數(shù)據(jù)中刪除，為此可刪除所找列的無關項ifexist('c:\matlabbk\sj.txt')loadc:\matlabbk\sj.txt[i,j]=size(sj)mu=mean(sj)sigma=std(sj)outliers=abs(sj-ones(i,1)*mu)>3*ones(i,1)*sigmasj(outliers)=[]elsesj=rand(40,3)end例3：輸入一個字符，若為大寫字母，則輸出其后繼字符，若為小寫字母，則輸出其前導字符，若為數(shù)字字符則輸出其對應的數(shù)值，若為其他字符則原樣輸出。

程序如下：

c=input('請輸入一個字符','s');ifc>='A'&c<='Z'disp(setstr(abs(c)+1));elseifc>='a'&c<='z'disp(setstr(abs(c)-1));elseifc>='0'&c<='9'disp(abs(c)-abs('0'));elsedisp(c);end例4：畫出下列分段函數(shù)所表示的曲面。[X,Y]=meshgrid(-3:0.1:3);T=X+Y;

ifT>1Z=0.54*exp(-0.75*X.^2-3.75*Y.^2-1.5*Y);

elseifT>-1andT<=1Z=0.7575*exp(-X.^2-6*Y.^2);

else

Z=0.5457*exp(-0.75*X.^2-3.75*Y.^2+1.5*Y);

endsurf(X,Y,Z)%plot3(X,Y,Z)%mesh(X,Y,Z)4、switch語句格式為：switch表達式case表達式1語句組1case表達式2語句組2……case表達式m語句組motherwise語句組m+1end例5:某商場對顧客所購買的商品實行打折銷售，標準如下(商品價格用price來表示)：price<200沒有折扣200≤price<5003%折扣500≤price<10005%折扣1000≤price<25008%折扣2500≤price<500010%折扣5000≤price14%折扣輸入所售商品的價格，求其實際銷售價格。price=input('請輸入商品價格');switchfix(price/100)case{0,1}%價格小于200rate=0;case{2,3,4}%價格大于等于200但小于500rate=3/100;casenum2cell(5:9)%價格大于等于500但小于1000rate=5/100;casenum2cell(10:24)%價格大于等于1000但小于2500rate=8/100;casenum2cell(25:49)%價格大于等于2500但小于5000rate=10/100;otherwise%價格大于等于5000rate=14/100;endprice=price*(1-rate)%輸出商品實際銷售價格try語句（試探性執(zhí)行語句）語句格式為： try語句組1catch語句組2endtry語句先試探性執(zhí)行語句組1，如果語句組1在執(zhí)行過程中出現(xiàn)錯誤，則將錯誤信息賦給保留的lasterr變量，并轉去執(zhí)行語句組2。這種試探性執(zhí)行語句是其他高級語言所沒有的。例6：矩陣乘法運算要求兩矩陣的維數(shù)相容，否則會出錯。先求兩矩陣的乘積，若出錯，則自動轉去求兩矩陣的點乘。程序如下：A=[1,2,3;4,5,6];B=[7,8,9;10,11,12];tryC=A*B;catchC=A.*B;endClasterr%顯示出錯原因3、循環(huán)結構（1）fort=e1:e2:e3

end

首先計算三個表達式的值，再將表達式1的值賦給循環(huán)變量，如果此時循環(huán)變量的值介于表達式1和表達式3的值之間，則執(zhí)行循環(huán)體語句，否則結束循環(huán)的執(zhí)行。執(zhí)行完一次循環(huán)之后，循環(huán)變量自增一個表達式2的值，然后再判斷循環(huán)變量的值是否介于表達式1和表達式3之間，如果滿足仍然執(zhí)行循環(huán)體，直至不滿足為止。這時將結束for語句的執(zhí)行，而繼續(xù)執(zhí)行for語句后面的語句。例1:已知求y的表達式，y=1/1+1/2^2+1/3^2+……＋1/n^2當n=100時，求y的值。程序如下：y=0;n=100;fori=1:ny=y+1/i/i;endy在實際MATLAB編程中，為提高程序的執(zhí)行速度，常用向量運算來代替循環(huán)操作：n=100;i=1:n;f=1./i.^2;y=sum(f)例如2：求定積分。程序如下：a=0;b=3*pi;n=1000;h=(b-a)/n;x=a:h:b;f=exp(-0.5*x).*sin(x+pi/6);fori=1:ns(i)=(f(i)+f(i+1))*h/2;ends=sum(s)事實上，MATLAB提供了有關數(shù)值積分的標準函數(shù)，實際應用中可直接調(diào)用這些函數(shù)求數(shù)值積分。

例3：fori=1:2x=input('x:')ifx>0y=x.^2+1subplot(1,2,1)plot(x,y)holdonelsey=x.^2-1subplot(1,2,2)plot(x,y)holdonendpauseend例4：用循環(huán)語句編寫M文件計算ex的值，其中x,n為輸入變量，ex的近似表達式為functiony=e(x,n)y=1;s=1;fori=1:ns=s*i;y=y+x^i/s;endyy=e(1,100)y=2.7183

for語句更一般的格式：for循環(huán)變量=矩陣表達式循環(huán)體語句end執(zhí)行過程是依次將矩陣的各列元素賦給循環(huán)變量，然后執(zhí)行循環(huán)體語句，直至各列元素處理完畢。實際上，“表達式1:表達式2:表達式3”是一個僅為一行的矩陣(行向量)，因而列向量是單個數(shù)據(jù)。例5：已知5個學生4門功課的成績，求每名學生的總成績。程序如下：s=0;a=[65,76,56,78;98,83,74,85;76,67,78,79;98,58,42,73;67,89,76,87];fork=as=s+k;enddisp(s');（2）While表達式<語句體>end其執(zhí)行過程為：若條件成立，則執(zhí)行循環(huán)體語句，執(zhí)行后再判斷條件是否成立，如果不成立則跳出循環(huán)例6：根據(jù)矩陣指數(shù)的冪級數(shù)展開式求矩陣指數(shù)。程序如下：X=input('X=');E=zeros(size(X));F=eye(size(X));n=1;whilenorm(F,1)>0E=E+F;F=F*X/n;n=n+1;endEexpm(X)%調(diào)用MATLAB矩陣指數(shù)函數(shù)求矩陣指數(shù)

循環(huán)的嵌套如果一個循環(huán)結構的循環(huán)體又包括一個循環(huán)結構，就稱為循環(huán)的嵌套，或稱為多重循環(huán)結構?？梢园凑涨短讓訑?shù)，分別叫做二重循環(huán)、三重循環(huán)等。處于內(nèi)部的循環(huán)叫作內(nèi)循環(huán)，處于外部的循環(huán)叫作外循環(huán)。在設計多重循環(huán)時，要特別注意內(nèi)、外循環(huán)之間的關系，以及各語句放置的位置，不要搞錯。

例7

用篩選法求某自然數(shù)范圍內(nèi)的全部素數(shù)。程序如下：m=input('m=');p=2:m;fori=2:sqrt(m)n=find(rem(p,i)==0&p~=i);p(n)=[];endp綜合實例用MATLAB求解問題時，一般要經(jīng)歷建模和編程兩個過程，只有在建模正確的前提下，方能得出正確的結果。

單自由度系統(tǒng)有阻尼自由振動1.建立計算模型

由動力學可知，單自由度有阻尼自由振動的振動方程為：無量剛化后有：其中

,上述方程的解為：其中

x0

表示初始位置，

v0

表示初始速度。

參數(shù)ωn＝10，x0＝1，v0＝0，計算的終止時間t=2。試求ξ從0.1到1運動方程的解，并畫出波形。2.MATLAB編程

編寫M文件ex1.m%首先清空MATLAB的工作空間clear;%給定初值

wn=10;tf=2;x0=1;v0=0;%計算不同的ξ值所對應的振型

for

j=1:10;eta(j)=0.1*j;wd(j)=wn*sqrt(1-eta(j)^2);%求振幅A

a=sqrt((wn*x0*eta(j)+v0)^2+(x0*wd(j))^2)/wd(j);綜合實例%求相位角phi=atan2(wd(j)*x0,v0+eta(j)*wn*x0);%設定自變量數(shù)組t

t=0:tf/1000:tf;%求過渡過程

x(j,:)=a*exp(-eta(j)*wn*t).*sin(wd(j)*t+phi);end%在同一個圖形窗口中繪制不同的ξ值所對應的振型

plot(t,x(1,:),t,x(2,:),t,x(3,:),t,x(4,:),...t,x(5,:),t,x(6,:),t,x(7,:),t,x(8,:),...t,x(9,:),t,x(10,:))gridon%新建一個圖形窗口，繪制三維網(wǎng)格圖

figuremesh(x)綜合實例綜合實例如果改變初始條件令x0＝0，v0＝1，其運動曲線實際上就是系統(tǒng)的脈沖過渡函數(shù)。綜合實例2、氣體分子運動的麥克斯韋分布曲線通過本例說明如何用復雜的數(shù)學公式繪制曲線。利用氣體分子運動的麥克斯韋速度分布律，求氯分子運動的速度分布曲線，并討論溫度T及分子量mu對速度分布曲線的影響。(1)建立計算模型麥克斯韋速度分布律為:其中，m---分子質量,m=mu/NA,mu---分子量,NA---阿伏加德羅數(shù)k---波爾茨曼常數(shù)T----氣體的絕對溫度v----分子速度綜合實例為研究單個參數(shù)的影響，先把麥克斯韋分布律編為一個函數(shù)子程序，以便重復調(diào)用，同時將常數(shù)項也放在子程序中。需要強調(diào)的是：子程序不得與主程序放在同一個M文件中，只能將子程序單獨做成M文件，并放在與主程序同一個工作路徑中。2.MATLAB編程

首先建立計算麥克斯韋分布律的子程序

mxw.mfunctionf=mxw(T,mu,v)%Thesubfunction

mxw.mofex2利用麥克斯韋速度分布律求分子的速度分布曲線的子程序%mu、v、T分別是分子量、分子速度和氣體的絕對溫度

k=1.381*10^(-23);%波爾茨曼常數(shù)

NA=6.022*10^23;%阿伏加德羅數(shù)m=mu/NA%分子質量

f=4*pi*((m/2*pi*k*T)).^(3/2).*v.*v.*exp(-m*v.^2./(2*k*T));綜合實例編寫主程序ex2.m

T=300;mu=28e-3;%給出T和mu的值

v=0:1500;%調(diào)出自變量數(shù)組y=mxw(T,mu,v);%調(diào)用子程序

plot(v,y,'r')%繪制分布曲線

holdon%為了看出不同的T和mu對曲線形狀的影響，再次給定T和mu，在同一幅圖中繪制分布律曲線的圖形

T=200;mu=28e-3;y=mxw(T,mu,v);plot(v,y,'b')holdonT=300;mu=2e-3;y=mxw(T,mu,v);plot(v,y,'g'))綜合實例3、方波的分解在連續(xù)信號系統(tǒng)中，方波可以用相應頻率的基波及其奇次諧波合成，這也是將方波展開為正弦級數(shù)的出發(fā)點。本節(jié)將演示這一現(xiàn)象。（1）建立計算模型

一個以原點為奇對稱中心的方波y(t)可以用奇次正弦波的疊加來逼近：

方波的寬度為π，周期為2π。

2.MATLAB編程建立M文件ex3.m%演示基波和奇次諧波合成方波t=0:0.1:10;%首先設定一個有101個點的時間數(shù)組%繪制頻率w＝1（f=1/2π）的正弦基波，并設置暫停

y=sin(t);plot(t,y)pause%疊加3次諧波，繪圖并設置暫停y=sin(t)+sin(3*t)/3;plot(t,y)pause%疊加1、3、5、7、9次諧波，繪圖并設置暫停y=sin(t)+sin(3*t)/3+sin(5*t)/5+sin(7*t)/7+sin(9*t)/9;plot(t,y)pause綜合實例%為了繪制三維曲面，需要將各次波形數(shù)據(jù)存儲為一個三維數(shù)組，因此需要重新定義y，重新編程，本例將求至19次諧波t=0:0.031:3.14;y=zeros(10,max(size(t)));x=zeros(size(t));fork=1:2:19x=x+sin(k*t)/k;y((k+1)/2,:)=x;end

pause%將各個波形疊合繪出，并設置暫停

plot(t,y(1,:),t,y(2,:),t,y(3,:),t,y(4,:),t,y(5,:),...t,y(7,:),t,y(8,:),t,y(9,:))pause%將各個波形繪制成三維網(wǎng)格圖mesh(y)pause綜合實例４、程序流控制交互包括兩個方面，輸入數(shù)據(jù)和暫停。

input命令格式：n=input（‘顯示信息字符串’）通過鍵盤輸入的數(shù)據(jù)可以是數(shù)值也可以是字符串。

pause命令：系統(tǒng)暫停狀態(tài)，也進入鍵盤主控狀態(tài)

keyboard進入鍵盤主控狀態(tài)，直接修改或輸入變量，用return可退出鍵盤控制狀態(tài)，系統(tǒng)繼續(xù)執(zhí)行后續(xù)操作．

echoechoon/off顯示/不顯示命令文件中的控制echo兩種狀態(tài)下切換

5.4函數(shù)變量及變量的作用域

控制輸入變量的函數(shù)（1）nargin：控制輸入變量的個數(shù)，使得在編程過程中，可以實現(xiàn)不定參數(shù)輸入的操作（2）varargin：實現(xiàn)不定數(shù)目輸入變量的函數(shù)設計。對函數(shù)的一切輸入變量均存儲在varargin命名的細胞數(shù)組中。例：當輸入一個矩陣時，求行列式值，當輸入二個矩陣時，求矩陣和Functionc=test69(a,b)if(nargin==1)c=det(a)elseif(nargin==2)c=a+bend5.4函數(shù)變量及變量的作用域Function[mathavg,engavg,chineseavg]=test610(varargin)l=length(varargin)mathsum=0engsum=0chinesesum=0fori=1:lmathsum=mathsum+varargin{i}(1)engsum=engsum+varargin{i}(2)chinesesum=chinesesum+varargin{i}(3)endmathavg=mathsum/lengavg=engsum/lchineseavg=chinesesum/l調(diào)用：test610（[60,70,80],[45,80,90]）控制輸出變量的函數(shù)Nargout：控制輸出變量的個數(shù)，使得在編程過程中，可以實現(xiàn)不定參數(shù)輸入的操作。Varargout：現(xiàn)不定數(shù)目輸出變量的函數(shù)設計。5.4函數(shù)變量及變量的作用域全局變量與局部變量

函數(shù)M文件的內(nèi)部變量是局部的，其他函數(shù)文件無法使用。為解決函數(shù)之間的傳遞數(shù)據(jù)的問題，就需要在函數(shù)文件執(zhí)行之前，定義變量為全局變量，其作用域是整個MATLAB的工作區(qū)，即全程有效。命令的格式：globalxyz舉例：z=α(x-1)2+β(y+1)2定義過程：函數(shù)M文件的內(nèi)容：functio