5.1MATLAB程序結構及其運行機制

5.2MATLAB程序的流程控制

5.3編寫MATLAB函數

5.4編程案例

第5單元MATLAB程序設計

MATLAB程序由主程序和函數(亦稱子程序)兩部分組成，如圖5-1所示。編程采用模塊化程序設計方法，用戶任務分解為多個子任務，子任務編寫成用戶自行命名的MATLAB函數。主程序設計成若干個模塊，每個模塊實現特定的數據處理目標，模塊調用MATLAB函數完成子任務，多個模塊順序銜接構成主程序，在CommandWindow中鍵入主程序的存盤文件名則可執行和完成全部任務。5.1MATLAB程序結構及其運行機制圖5-1MATLAB程序的結構與數據流5.2.1for循環

(1)利用循環向量化技術設計for循環，計算數列｛135…｝的前15項累加和與累乘積。程序如下：

5.2MATLAB程序的流程控制clc;closeall;clearall;

x=0;y=1; %賦初值

fork=1:2:2*15-1; %創建循環變量

x=x+k;y=y*k; %計算累加和與累乘積

end;

formatshorte;Results=[xy] %計算結果

程序執行結果如下：

Results=

2.2500e+0026.1903e+015

(2)利用陣列預分配技術設計for循環，計算下面離散系統在正弦輸入下的響應yk。

程序如下:

clc;closeall;clearall;tic %初始化和啟動秒表計時器

space=0.001;t=0:space:16;u=sin(2*pi*t); %計算正弦輸入

n=length(u); %計算u的數據個數

y=zeros(n,1); %按u的數據個數n、實施陣列y的預分配

y(1)=2*u(1);y(2)=0.75*y(1)+2*u(2);%計算前2個不能遞推計算的響應

fork=3:n; %創建循環變量

y(k)=0.75*y(k-1)-0.125*y(k-2)+2*u(k); %用遞推公式計算響應y

end;

Time=toc %測定程序運行時間

plot(t,y) %響應y的計算結果展示

程序執行結果如下，圖形如圖5-2所示。

Time=

5.1431e-003圖5-2離散系統的時間響應y(t)5.2.2while循環

while循環程序設計應注意兩點，一是循環變量累加，二是循環變量關系表達式，兩者構成循環控制邏輯。

用關系表達式設計while循環，計算函數在區間上的面積。程序如下：

clc;closeall;clearall;tic %初始化和啟動秒表計時器

space=pi/30;theta=0:space:pi;y=sin(theta); %生成被積函數數據

x=0;k=1; %給面積變量x和循環變量k賦初值

whilek<=length(y); %用下標k作循環變量，設計關系表達式

%whilek<length(y)+1; %可選的循環變量關系表達式x=x+y(k)*space; %累計計算面積和

k=k+1; %循環變量增1累加

end;

time＝toc

disp('k='),disp(k)

disp('aera='),disp(x)

程序執行結果：

time=

0.0065

k=

32

aera=

1.99825.2.3if-else-end結構

if-else-end結構關系表達式根據檢驗結果有條件地執行程序。

(1)編寫一個MATLAB函數用于判斷一個小于50的數是否為素數，并存盤為prime.m。程序如下：

functiony=prime(x)

ifrem(x,2)==0&x>2;

y='不是素數,';

elseifrem(x,3)==0&x>3;

y='不是素數,';

elseifrem(x,5)==0&x>5;

y='不是素數,';

elseifrem(x,7)==0&x>7;

y='不是素數,';

else;

y='是素數，';

end;編寫主程序指定要判斷的多個數字(向量x)，存盤名為myprime01.m。程序如下：

clc;closeall;clearall;

x=[123456789101112131415]; %指定要判別的數字

n=length(x);p=[]; %判別結果變量p賦初值

fori=1:n; %設計循環逐個判別

p=[p[num2str(x(i))prime(x(i))]]; %判別結果放到一個字符串變量p里

end;

disp(p);%顯示判別結果

向量素數判斷主程序的執行結果如下：

1是素數,2是素數,3是素數,4不是素數,5是素數,6不是素數,7是素數,8不是素數,9不是素數,10不是素數,11是素數,12不是素數,13是素數,14不是素數,15不是素數,

(2)將判斷一個矩陣各個元素是否為素數的任務編寫成一個MATLAB函數，并存盤為primematrix.m。程序如下：

functiony=primematrix(x)

x1=rem(x,2)==0&x>2;

x2=rem(x,3)==0&x>3;

x3=rem(x,5)==0&x>5;

x4=rem(x,7)==0&x>7;

y=find((x1+x2+x3+x4)==0);

編寫主程序指定要判斷的多個數字(矩陣x)，存盤名為myprime02.m。程序如下：

clc;closeall;clearall;

x=[12345;678910;1112131415]%指定要判別的矩陣

p=[primematrix(x)]'%判別結果為向量p，其元素值為矩陣x中素數元素的序號矩陣素數判斷主程序的執行結果如下：

x=

12345

678910

1112131415

p=

134579135.2.4switch-case-end結構

switch-case-end結構根據狀態表達式的值執行不同語句塊，相當于多條if塊的嵌套使用。

將判斷一個數是否為素數的任務編寫成一個MATLAB函數并存盤為oneprime.m。程序如下：

functiony=oneprime(x)

u=(rem(x,2)==0&x>2)|(rem(x,3)==0&x>3)|(rem(x,5)==0&x>5)…

|(rem(x,7)==0&x>7);

switchu;

case0;

y='是素數,';

case1;

y='不是素數,';

otherwise;%這部分語句塊在本函數可省略

y='是素數,';

end;編寫主程序指定要判斷的多個數字(向量x)，存盤名為myprime03.m。程序如下：

clc;closeall;clearall;

x=[123456789101112131415];%指定要判別的數字

n=length(x);p=[];%判別結果變量p賦初值

fori=1:n;%設計循環逐個判別

p=[p[num2str(x(i))oneprime(x(i))]];%判別結果放到一個字符串變量p里

end;

disp(p);%顯示判別結果

主程序執行結果如下：

1是素數,2是素數,3是素數,4不是素數,5是素數,6不是素數,7是素數,8不是素數,9不是素數,10不是素數,11是素數,12不是素數,13是素數,14不是素數,15不是素數,5.3.1MATLAB函數的結構及其調用機制

MATLAB函數是一種特定結構的程序，在MATLAB環境里的程序與一般程序一樣都是可視的文本，因其存盤名的擴展名為“.m”，亦稱作M-file函數。

與一般程序不同，M-file函數不能在CommandWindow中通過簡單鍵入它的存盤基本名執行，而是采用按M-file函數定義行所指定的函數格式鍵入函數表達式的方法調用，主程序調用函數亦應遵循指定的函數格式。M-file函數的標準結構如圖5-3所示。5.3編寫MATLAB函數圖5-3M-file函數的結構在Editor窗口編寫M-file函數，存盤名由“基本名+擴展名”構成，存盤時僅僅鍵入基本名，擴展名默認為“.m”，函數調用時使用存盤的基本名和指定的函數格式。注意，MATLAB7.1版本以后存盤函數基本名與函數定義行指定的函數名可以不一致，而早先的版本要求存盤基本名與函數定義行的指定函數名必須一致。下面的函數示例計算(y,outp)=f(x,inp)，其中輸入參數x為自變量，輸入參數inp為指定擬合的多項式階數，輸出參數y為因變量，且

，輸出參數outp為依據(x,y)擬合的inp階多項式。程序如下：

function[y,outp]=myfun(x,inp);

y=100*sin(x).*exp(-x.^2);

outp=polyfit(x,y,inp);命名所編寫的M-file函數并存盤。(如存盤為myfun.m)

在CommandWindow窗口鍵入函數表達式，且預先為輸入參數賦值；在主程序里編寫函數表達式調用函數，且預先為輸入參數賦值。命令如下：

x=1:7;order=3;

[y,p]=myfun(x,order)

執行函數調用的結果如下：

y=

30.95601.66540.0017-0.0000-0.0000-0.00000.0000

p=

-0.813611.6055-51.532969.80295.3.2編寫匿名函數

匿名函數指在主程序中以表達式方式定義的函數，下面示例展示它的格式和調用方法，程序如下：

quade=@(x,a,b,c)a*x.^2+b*x+c; %定義匿名函數并賦給變量quade

x=[123];a=2;b=-1;c=5; %輸入參數賦值

y=quade(x,a,b,c)%以調用變量quade調用匿名函數，注意輸入參數的位置和意義

主程序執行結果如下：

y=

611205.3.3編寫標準M-file函數

標準M-file函數，指函數體中的語句不調用其它任何函數，均采用直接編寫的算式。

下面函數在x的多個點上計算二次函數的值，存盤名為stdfun.m。程序如下：

functiony=stdfun(x,a,b,c)

y=a*x.^2+b*x+c;

編寫調用函數stdfun的主程序。程序如下：

clc;closeall;clearall;%初始化

x=[123];a=2;b=-1;c=5;

y=stdfun(x,a,b,c)

主程序調用函數的執行結果如下：

y=

611205.3.4編寫私人函數

私人函數是一種特殊的M-file函數，特別放置在名為private的文件夾里，僅僅能被父文件夾(上一級文件夾)里的函數調用。父文件夾之外的函數無法調用私人函數，故私人函數可以與其他文件夾的函數名相同，且一個函數調用M-file函數時優先搜尋私人函數。若用戶想限制某個M-file函數的使用權限或隱蔽調用，可將其放置在private文件夾里。

假定用戶文件夾mydir在MATLAB的搜索路徑上，mydir的一個子文件夾命名為private，則private里的函數(稱作私人函數)僅僅能被mydir里的函數調用。注意，是父文件夾里的函數調用私人函數，而不是父文件夾里的主程序調用。創建用戶自己的文件夾mydir，編寫下面函數myprog.m并存放在mydir文件夾下，該函數設計為調用私人函數sin，操作File菜單的setpath項將mydir文件夾添加到MATLAB搜索路徑。

functiony=myprog(x)

y=sin(x)

在mydir文件夾內創建private文件夾，編寫私人函數sin.m并存放在private文件夾里。程序如下：

functiony=sin(x);

y=x.^2;%故意編寫一個錯誤表達式編寫下面主程序mytest.m并可存放在任何MATLAB搜索路徑的文件夾里。程序如下：

x=[123;456]

y_common=sin(x) %此表達式將調用MATLAB系統的函數sin.m

y_private=myprog(x) %此表達式將調用mydir文件夾下函數myprog.m

%myprog.m函數將調用private文件夾下私人函數sin.m

試比較如下主程序的執行結果：

x=

123

456

y_common=

0.8415

0.90930.1411

-0.7568-0.9589-0.2794

y_private=

149

1625365.3.5編寫形如“主函數+子函數”的M-file函數

如果一個M-file函數的函數體內調用一個或多個其它M-file函數，則稱其為主函數，被調用的M-file函數稱作子函數。在用戶任務里，主函數被主程序或命令行調用，而子函數僅僅被主函數或其它子函數調用。

下面的主函數rank調用svd、nargin、max、size、sum等5個子函數，以允差tol計算矩陣A的秩，主函數rank的存盤名為rank.m。程序如下：functionr=rank(A,tol)%RANKMatrixrank.s=svd(A);ifnargin==1;tol=max(size(A)')*max(s)*eps;end;r=sum(s>tol);編寫主程序計算矩陣x的秩。程序如下：

x=[1524;3679;5293];

rx=rank(x,1e-005)

主程序調用主函數rank的執行結果如下：

rx=

35.3.6編寫形如“主函數+嵌套函數”的M-file函數

“主函數+嵌套函數”是一種特殊結構的M-file函數，主函數供主程序或其它函數調用，嵌套函數按包含關系逐層嵌套在主函數內。一個嵌套函數可調用任何外部的M-file函數，但它自身只能在主函數內部被調用，且遵循下述規則：(1)主函數或父嵌套函數可調用子嵌套函數；(2)主函數或嵌套函數不能越層調用深層嵌套函數；(3)底層嵌套函數可調用高層嵌套函數，包括越層調用；(4)一個嵌套函數可調用同一層其它嵌套函數。

編寫下面嵌套函數，存盤為mynested.m。

function[y,p]=mynested(x,a,b,c,order)

y=computing(x,a,b,c);

[p,x1,y1]=polyfiting(x,y,order);

ploting(x,y,x1,y1);

functiony_value=computing(x,a,b,c)

y_value=a*x.^2+b*x+c;

end

function[p_poly,x1,y1]=polyfiting(x,y,order)

p_poly=polyfit(x,y,order);

[x1,y1]=interposing(x,y);

function[x1,y1]=interposing(x,y)

x1=min(x):0.01:max(x);

y1=interp1(x,y,x1);

end

end

functionploting(x,y,x1,y1)

plot(x1,y1,'--',x,y,'o');

end

end

編寫下面主程序調用嵌套函數。

clc;closeall;clearall;

x=linspace(-1.5,3,21);a=2;b=-1;c=5;order=3;

[y,p]=mynested(x,a,b,c,order)主程序調用嵌套函數結果之一如下：

y=

Columns1through11

11.00009.52628.25507.18626.32005.65635.19504.93634.88005.02635.3750

Columns12through21

5.92636.68007.63628.795010.156311.720013.486315.455017.626320.0000

p=

0.00002.0000-1.00005.0000

主函數調用嵌套函數結果之二如圖5-4所示。圖5-4主程序調用嵌套函數的結果之二5.3.7編寫函數的函數

函數的函數(functionfunctions)是一種供求根(Zerofinding)、優化(Optimization)、積分(Quadrature)和解常微分方程(Ordinarydifferentialequations)等工具函數調用的函數。詳見第4單元積分和解常微分方程的函數編程，亦可參見第8單元優化設計的函數編程。5.3.8全局變量和局部變量

在M-file函數中定義一些變量x1,x2,…為全局變量，使用語句globalx1x2…。如果M-file函數中的變量沒有定義為全局變量(GlobalVariable)，則默認全是局部變量(LocalVariable)。外部程序或外部M-file函數不能直接調用一個函數內的局部變量，故變量名可以與外部變量同名，但卻可直接調用全局變量。編程采用全局變量，遵循先定義后使用的規則。盡量避免采用全局變量。

編寫含全局變量a、b、c的M-file函數quading.m。程序如下：

functiony_value=quading(x)

globalabc;

y_value=a*x.^2+b*x+c;

end

編寫含同名全局變量a、b、c的主程序，實現主程序與函數之間的參數傳遞，程序如下：

clc;closeall;clearall;

globalabc;

x=linspace(-1.5,3,21);a=2;b=-1;c=5;order=3;

y=quading(x)主程序調用函數的執行結果如下：

y=

Columns1through11

11.00009.52628.25507.18626.32005.65635.19504.93634.88005.02635.3750

Columns12through21

5.92636.68007.63628.795010.156311.720013.486315.455017.626320.00005.3.9函數句柄

無論是系統的M-file函數或是自定義的M-file函數，都可以用函數句柄(FunctionHandle)來代表，也可用一個自定義的變量名替代，程序調用函數被調用函數句柄代替，這樣可使程序變得簡潔和更容易理解。

創建函數句柄并運用它。程序如下：

clc;clearall;closeall;

sqr=@(x)x.^2; %為匿名函數f(x)=x2創建一個函數句柄sqr

x=0:pi/3:2*pi;

y=sqr(5) %調用函數句柄sqr計算函數f(x)=x2的值

sqr=@sin; %為系統的M-file函數sin創建一個函數句柄sqr

a=sqr(x) %調用函數句柄sqr計算函數f(x)=sin(x)的值

b=quad(sqr,0,pi) %調用函數句柄sqr計算函數f(x)=sin(x)在區間[0,π]上的積分上面程序的執行結果如下：

y=

25

a=

00.86600.86600.0000-0.8660-0.8660-0.0000

b=

2.0000利用函數句柄技術編寫一個“主函數+嵌套函數”結構的M-file函數，主函數命名并存盤為get_plot_handle.m，主函數的輸入參數為plot圖所需的線型設置字符串變量LS、線寬數值變量LW、點標記邊界顏色字符串變量MEC、點標記表面顏色字符串變量MFC和點標記尺寸數值變量MS。嵌套函數命名為draw_plot，它的輸入參數為x,y，函數體中的plot函數調用主函數變量LS、LW、MEC、MFC和MS。程序如下：

functionh=get_plot_handle(LS,LW,MEC,MFC,MS)

h=@draw_plot;

functiondraw_plot(x,y)

plot(x,y,LS,'LineWidth',LW,...

'MarkerEdgeColor',MEC,...

'MarkerFaceColor',MFC,...

'MarkerSize',MS)

end

end采用函數句柄技術編寫主程序調用上面的函數。程序如下：

clc;clearall;closeall;

h=get_plot_handle('--rs',2,'k','g',10);%創建函數句柄h加載圖形設置數據

x=-pi:pi/10:pi; %計算繪圖數據x

y=tan(sin(x))-sin(tan(x)); %計算繪圖數據y

h(x,y) %調用函數句柄h，它按加載的圖形設置數據繪圖

主程序執行的結果如圖5-5所示。圖5-5采用函數句柄技術編程的繪圖結果5.4.1分析圖像的灰度分布

編寫顯示原圖的M-file函數original_image.m。程序如下：

functionoriginal_image(image)

%ShowaOriginalImage

imshow(image)5.4編程案例編寫顯示圖像灰度直方圖的M-file函數gray_histogram.m。程序如下：

functiongray_histogram(image,N);

%ShowtheGrayscaleHistogramofaImage

ifnargin==1;

imhist(image);boxoff;

else;

imhist(image,N);boxoff;

end;

編寫主程序myimage_analysis.m，實現原圖和灰度頻數直方圖的顯示。程序如下：

clc;clearall;closeall;

char1='D:\users\pictures\';char2='sample01.jpg';

h1=@original_image;h2=@gray_histogram;

xx=imread([char1char2]);yy=rgb2gray(xx);

subplot(2,2,1),h1(xx);title('OriginalImage')

subplot(2,2,2),h1(yy);title('GrayImage')

subplot(2,2,3),h2(yy);title('GrayscaleHistogram')

subplot(2,2,4),h2(yy,99);title('GrayscaleHistogram')

主程序的執行結果如圖5-6所示。圖5-6案例主程序的執行結果5.4.2懸掛機組機具跟蹤性分析

輪式拖拉機與鏵式犁懸掛連接組成懸掛犁耕機組，當拖拉機轉彎時，機具末端外側偏離理想軌跡的程度稱作機具跟蹤性。機組轉彎參數滿足下面的懸掛跟蹤方程：

試求機具偏角與拖拉機前輪轉角的函數關系=f?(θ)(相關參數如圖5-7所示)。圖5-7懸掛犁耕機組和轉彎參數

解:懸掛跟蹤方程是一個隱函數方程，不能直接獲得的解析表達式。因此，擬采用數值解法求出?

與??的n對觀察值，再用多項式擬合求出的近似函數關系。

算法流程設計如圖5-8所示。圖5-8懸掛犁耕機組機具跟蹤性分析程序的算法流程編寫懸掛跟蹤方程求根的M-file函數myequation.m，即給定θ解出。程序如下：

function[phi,fval]=myequation(theta)

L=16*0.3048;h=10*0.3048;b=4*0.3048;

n=length(theta);tol=1e-5;

m=floor(pi/(2*tol));

phi0=linspace(0,pi/2,m);

phi=zeros(n,1);index=zeros(n,1);fval=zeros(n,1);

fori=1:n;

f0=((L+h)*cos(phi0)-b*sin(phi0))*sin(theta(i))-L*sin(theta(i)+phi0);

[fval(i),index]=min(abs(f0));

ph