10/30/2006

第八章 多項式之應用

對於多項式MATLAB也提供許多指令可供運算,相關的指令如下表:

函數名稱 說明
conv(p,q) 兩多項式相乘
[q,r]=deconv (num,den) 兩多項式相除,q為商,r為餘數,num分子,den為分母
roots 求多項式之根
poly(r) 將根轉為多項式.
polyval(p,x) 計算多項式之值,p為多項式之係數矩陣,x為變數值,可為矩陣型式
polyvalm(p,X) 計算多項式之值,p為多項式之係數矩陣,X為變數值,為方矩陣型式
residue 部份展開式之餘數.
polyfit(X,Y,n) 多項式資料回歸,(X,Y)為對應資料,n多項式最高次方
polyder (p) 多項式微分
polyint(p,K) 多項式之積分,K為常數


多項式中,主要以其係數組成一向量作為運算之基礎。多項式之通式可表示如下:


F(x)=a1xn+a2xn-1+...+an-1x2+anx+an+1


其中x 為變數,n為其最高之階,為變數x之次方。而a1, a2 ,…,an , an+1等則分別係數,可以利用向量矩陣表示。以上項之通式為例,其向量p為:


p=[a1, a2 ,…,an , an+1]


代表F(x)這個多項式。多項式間之運算則以其p向量為代表。例如,兩多項式之相加,設其向量p=[20 -7 5 10],另一為q=[1 8 1 -6],則相加後應為m=[21 1 6 4],亦即其多項式和為21X³+X²+6X+4。此處若兩向量之大小不同,則較少的前面元素應以零補齊。兩多項式相減時亦同。多項式若乘以一個常數,則直接以常數乘上每項係數或向量中之每個元素即可。但若為兩個非常數之多項數相乘或相除,則必須藉助MATLAB之指令conv或 deconv來運算才行。這兩個指令均有兩項輸入,分別為兩多項式之向量。例如兩多項式之向量分別為p,q:


>>p=[1 5 3 4 1] %代表多項式X+5X3+3X 2+4X+1
p = 1 5 3 4 1

>>q=[8 3 -5 6] %代表多項式8X3+3X2-5X+6
q = 8 3 -5 6

兩者相乘後,設為M,則M應為:

多項式乘法conv



>>M=conv(p,q)
M =
8 43 34 22 35 1 19 6


此多項式為

8X7+43X6+34X5+22X4+35X3+X2+19X+6。


多項式除法deconv


今若將M除以q,利用deconv函數指令,其結果應為:


>>[s,r]=deconv(M,q)
s =
1 5 3 4 1
r =
0 0 0 0 0 0 0 0

此s向量與p向量相同,r則為除後之餘數,目前為零,因為完全除盡。例如:

>>[s,r]=deconv(p,q)
s =
0.1250 0.5781
r =
0 0 1.8906 6.1406 -2.4688

此時之商s為0.1250X + 0.5781;其餘項為

1.8906X3+ 6.1406X2 -2.4688


上述多項式相乘除之指令conv與deconv並不需要相同的向量長度。即使除法,其前項為分子,後項為分母,分子之長度通常比分母長,但並不必要如加法或減法,需在前面加零元素,以湊成相同的長度。

多項式求值polyval


其次為多項式值之計算可用polyval(p,x)指令計算多項式之值,其中p為多項式之係數向量,x為變數值,可為向量矩陣型式:


>>x=0:5,p=[1 -3 3 6]
x =
0 1 2 3 4 5
p =
1 -3 3 6

>>f=polyval(p,x)
f =
6 7 8 15 34 71


上式中p代表多項式f(x)之向量,其所得值為f(0)、f(1),…f(5)分別為f向量之元素值。運算時,亦可合併為一個指令執行,只要位置放對就行:


>>f=polyval([1 -3 3 6],[0:5])
f =
6 7 8 15 34 71

8.1多項式之根

令一多項式等於零,其x值即為其根。具有n階的多項式,應存在有n個根,但這些根可能為實數,也可能為複數或重根。若根為複數,只要其向量元素為實數,則其複數根常成對存在,或稱為共軛根,其型式為a±bi。

roots


求根的指令為roots(p),p為多項式之向量矩陣。如:


>>roots([2 -2 -1 0 1])
ans =
1.0000 + 0.0000i
1.0000 - 0.0000i
-0.5000 + 0.5000i
-0.5000 - 0.5000i

結果得到兩個相同的實根1及-5±0.5i。在多項式中,最簡單的二次式ax²+bx+c=0之根可利用公式得解,其型式如下:

x={(-b ±(b²-4ac)(1/2)}/2a

若二次之型式為ax²+2bx+c=0,則其根之型式可改變為:

x={(-b ±(b²-ac)(1/2)}/a

設a=5,b=10, 則上式根可以計算如下:


>>a=2;b=10;c=12;
>>x1=(-b+sqrt(b.*b-4*a.*c))/(2*a)
x1 = -2

>>x2=(-b-sqrt(b.*b-4*a.*c))/(2*a)
x2 = -3

若用roots指令求解,則


>>roots([2 10 12])
ans =
-3.0000
-2.0000

所得之答案相同。

poly


有根時,也可以利poly(r)指令將其轉回多項式之型式。以前面所得之兩個實根1及虛根 -5±0.5i為例,其對應多項式應為:


>>poly([1, 1, 0.5+0.5i, -0.5-0.5i])
ans =
1.0000 -1.0000 -0.5000 0 0.5000

其結果若分別乘以2後,應與前例之設定相同。由此可知,用根反求多項式係數時,其結果可能為整數倍數,其比例應一致的。而且,在poly(r)之輸入項r可為行向量或列向量,其結果應相同。

8.2 非線性函數求根法(1)

除多項式之尋根外,有些函數屬於非線性型式,其尋根法自有不同。在matlab中,有fzero、fminbnd及 function_handle可以運用。茲介紹如下:

FZERO 尋找函數根之指令


fzero之語法如下:

 x = fzero(fun,x0)
 x = fzero(fun,x0,options)
 [x,fval,exitflag,output] = fzero(fun,x0,options)


通常要對任一個函數求其等於零之根,有如大海撈針,尤其函數若在零點時存在有許多根時,其所找出之根有時並不符合所望。為此,必須先給附近值或概略區間,使程式能順利找到所要的根。fzero指令之功能就是利用x0為第一個猜測值進行求根。若x0為常數,程式會先找到一個具有正負值區間但包含x0值之根。若設定根之範圍,則可改為矩陣之型式,不過所設定之範圍若其對應函數值無法產生正負符號,程式可能會給一個錯誤的信息。

在輸出參數中,除得到根x之外,其對應之函數值照理為零,但若有差異,亦可利用第二個輸出參數fval。若在設定之範圍內得到NaN或Inf時,表示無解,但也可能以複數為最後的答案。FUN之參數則代表函數名稱,包括自訂函數在內。其名稱之前則需加@之符號。例如:

>> X = fzero(@sin,3)
X =
3.1416

>>> x = fzero(@cos,[1 2])
x =
1.5708

函數也可以使用匿名函數,例如,若有一個匿名函數為 @(x) sin(3*x/5):

>> X = fzero(@(x) sin(3*x/5),5)
X =
5.2360

函數也可使用自訂函數,例如下面之自訂函數demofun:

function f=demofun(x)
f = x.^3-100;

>> x=fzero(@demofun,1)
x =
4.6416


若自訂函數有兩個以上輸入參數,則可採用匿名函數之型式,但只能針對其中一個當變數,其餘需事先設定為常數,例如:


>> x=fzero(@(x) demofun(x),1)
x =
4.6416


使用匿名函數時,其函數名稱也可使用函數握把,但不必加上@,例如:

>>f = @(x)x.^3-5*x-6;

>> z=fzero(f,3)
z =
2.6891


參數options則以optimset函數產生,主要是設定尋根的過程參數,其中包括Display, TolX, FunValCheck及OutputFcn等項目。此外,exitflag輸出參數則是執行結果之代碼,其中:





1 FZERO找到適當值
-1運算依output函數之設定停止。
-3在設定之區間,其結果為NaN或Inf。
-4找到結果為複數型式。
-5 FZERO可能接近單質點(singular point)

8.3多項式分母之展開式

多項式餘式residue


多項式分母若無重根的情況,應可進行因式分解,並求得其係數,此可以利用residue指令求解,其語法如下:

[r,p,k] = residue(b,a)
[b,a] = residue(r,p,k)


其中,b與a分別為兩個多項式。其中a為一分式之分母,b該分式之分子,兩者均為向量型式;而[r,p,k]等則為行向量,分別代表展開後之分子、分母及餘數之係數,若能除盡則k項應為零。r與p分別為餘數與極數,其個數應相等,但應比a之個數少一。

當極點(根)均相異的狀況下,值若p中有相同值時,其解的型式必須稍作變化,例如:

b =[ 3 -6 4]; a =[ 1 -5 8 -4];
>> [r,p,k]=residue(b,a)
r = 2
4
1
p = 2
2
1
k = []

其結果應解釋為如下之型式:

b(x)/a(x)=(3x²-6x+4)/(x³-5x²+8x-4)=2/(x-2)+4/(x-2)²+1/(x-1)

RESIDUE這個指令也可以採用逆向方式反求A與B之多項數,只要將其參數反向即可,如:

[B,A] = RESIDUE(R,P,K)

其中RPK之定義如前述,但此時成為輸入項,輸出為A與B,但應注意其對應位置。下面為上述之反向例:

>> [b,a]=residue([2 4 1],[2 2 1],[])
b = 3 -6 4
a = 1 -5 8 -4

其所得結果與前述相同。

8.4 非線性函數求根法-fminbnd

FMINBND 非線性函數最小值


fminrnd之功能是在特定區間尋找單一變數之函數之最小值。其指令型式如下:

x = fminbnd(fun,x1,x2)
x = fminbnd(fun,x1,x2,options)
[x,fval,exitflag,output] = fminbnd(...)



上述指令fminbnd之功能在於針對函數fun,在設定的範圍x1<x<x2中找尋其最小值。其值置於左邊之x中,而函數之對應值則置於fval內。其中fun為函數之握把,或函數之名稱。若為函數名稱則須在名稱前加@符號。


>> X = fminbnd(@sin,3,4)

X =

3.9999


也可以使用匿名函數:

>> f = @(x)x.^3-2*x-5;
>> [x,f] = fminbnd(f, 0, 2)
x =
0.8165
f =
-6.0887

>> x = fminbnd(@(x) sin(x)+3,2,5)

x =

4.7124



在指令中之第四參數項options之參數則可利用optimset函數依下列參數進行設定:








Display顯示之層次:'off'不顯示;'iter'顯示每次回圈之結果; 'final' 只顯示最後結果;'notify' (預設) 若函數不收斂時顯示。
FunValCheck檢查目的函數是否正確,'on' 目標函數之結果為複數或 NaN時顯出警告;'off' 不顯示警告。
MaxFunEvals函數容許之最大值。
MaxIter容許最大迴圈數。
OutputFcn指出使用者設定之呼叫定義函數。
TolX x值之容許度。


上述參數之設定可利用optimset函數,例如:

>> [X,FVAL,EXITFLAG] = fminbnd(@cos,3,4,optimset('TolX',1e-12,'Display','off'))

X =
3.1416
FVAL =
-1
EXITFLAG =
1


輸出之第三參數為運算之狀態,與fzero中之第三參數相同:






1fminbnd指令依據TolX收斂至x。
0最大迴圈數。
-1由 output函數決定終止。
-2邊界不符合 (x1 > x2)。


輸出之第四項為OUTPUT,其內容如下:






output.algorithm使用之運算法
output.funcCount函數計算之次數
output.iterations迴圈之次數
output.message離開信息

8.5多項式曲線適配法

多項式曲線適配法可以利用polyfit完成。利用現有之資料組合[x,y]以求最佳之多項式曲線,其語法如下:


p = polyfit(x,y,n)
[p,S] = polyfit(x,y,n)
[p,S,mu] = polyfit(x,y,n)


其中[x, y]為資料組,n為多項式之最高階數,p則是所得多項式之係數向量。S則為結構矩陣,其下有R、df、normr等欄位,代表QR之分解結構、自由度及副變量。

範例:


這個例子是利用現有資料適配erf(x)函數,並利用x之多項式表示。這個例子雖然不很適當,因為erf(x)函數有其特定範圍,一般多項式則沒有範圍。故利用資料適配的結果可能不會很理想。

開始時先產生x之向量點,該其均勻分佈於一個區間,然後使用erf(x)函數進行估計:

x = (0: 0.1: 2.5)';
y = erf(x);

設多項式最高階n=6,則

p = polyfit(x,y,6)

得到之結果,多項式係數為:

p =

0.0084 -0.0983 0.4217 -0.7435 0.1471 1.1064 0.0004


為瞭解所得之多項式與實際值有多接近,可以使用polyval這個求值的函數,同時將其作成表格作比較,最後一項則為誤差:


f = polyval(p,x);

table = [x y f y-f]

table =

0 0 0.0004 -0.0004
0.1000 0.1125 0.1119 0.0006
0.2000 0.2227 0.2223 0.0004
0.3000 0.3286 0.3287 -0.0001
0.4000 0.4284 0.4288 -0.0004
...
2.1000 0.9970 0.9969 0.0001
2.2000 0.9981 0.9982 -0.0001
2.3000 0.9989 0.9991 -0.0003
2.4000 0.9993 0.9995 -0.0002
2.5000 0.9996 0.9994 0.0002

在這個區間內,看起來其適配情形還不錯,可以到小數三、四位。但不要高興太早,因為超出此範圍後,多項式會如脫韁之馬,很快就露出馬腳。利用下面程式內容,可以立即用圖形進行比較其誤差之大小,結果則請讀者自行執行。


x = (0: 0.1: 5)';
y = erf(x);
f = polyval(p,x);
plot(x,y,'o',x,f,'-')
axis([0 5 0 2])

8.6 非線性函數求根法-fminsearch

FMINSEARCH指令


此為尋求函數之最小值位置之另一個指令,或稱為尼德-米(Nelder-Mead)法。它可應用於多維非線性函數,不必使用範圍界定。指令之格式如下:

x = fminsearch(fun,x0)
x = fminsearch(fun,x0,options)
[x,fval,exitflag,output] = fminsearch(fun,x0,options)


fminsearch指令起始值為x0,尋找其附近fun函數之最小值,結果置於x。x可為常數,向量或矩陣。同理,options參數可由optimset函數設定。其項目包括 Display, TolX, TolFun, MaxFunEvals, MaxIter, FunValCheck及OutputFcn,讀者可以參閱手冊。

>> X = fminsearch(@cos,3)

X =
3.1416

---------------------------
>> f=@(x) cos(x)+sin(x);

>> X = fminsearch(f,[5])

X =

3.9270

>> X = fminsearch(@(x) cos(x)+sin(x),[5])

X =

3.9270


8.7多項式矩陣值

前述之polyval可計算一個多項式值,若輸入值為矩陣型式,則可使用polyvalm求值,其語法如下:


Y= polyvalm(p, X)


其中,p為多項式之係數向量,X為其變數值,以矩陣型式表示,但必須為方矩陣。
例如:

>> p=[1 5 12 40 15]
p =
1 5 12 40 15
>> X=magic(4)
X =
16 2 3 13
5 11 10 8
9 7 6 12
4 14 15 1
>> m=polyval(p,X)
m =
89743 199 459 42109
1765 23203 16615 7759
11553 4999 3063 31599
943 55063 70815 73

>> m=polyvalm(p,X)
m =
395489 381218 383866 387530
382538 390345 389154 386066
386506 387834 388145 385618
383570 388706 386938 388889

10/28/2006

第七章 資料之輸入與輸出

資料之輸入及輸出同等重要。輸入與輸出均有特定的指令,有些應用於檔案,有些則直接由鍵盤輸入。在輸出的指令中,比較正式使用的為fprintf函數,其型式係沿用C語言,有些設定與該語言相通;比較簡單的型式則如disp指令,主要應用於螢幕之顯示。

輸出指令則以fscanf為主,其參數及格式與fprintf大體相同,而螢幕常用的指令則是input函數。詳細的分類可參考表5.1。

表7.1輸入輸出相關指令















輸出入指令說明
input 提供字串,使使用者輸入資料於一變數中
keyboard將執行權交至鍵盤
fscanf 按照特定的式自某輸出裝置輸出,包括輸出於螢幕上。
fread 讀入二位元資料
fprintf按照特定格式,印出資料
disp(x)輸出矩陣內容或文字內容,但不含變數名稱
menu自訂一個選擇菜單,經選項後決定其值
ginput 利用滑鼠之位置輸入座標值
return 結束返回
pause 暫停執行,按任意鍵繼續
fopen, fclose開啟檔案、關閉檔案
fwrite 寫入二位元資料

7.1 INPUT('message')

input之輸入指令主要應用於鍵盤輸入,然後按enter 確定輸入之內容。這些輸入之內容將設定給左邊的變數。


>>age=input('請問芳齡:')
請問芳齡: 60
age = 60


這個指令執行時,會將訊息打在螢幕上,然後等待使用者輸入任何資料。若沒輸入任何資料,僅按enter鍵,則age將會為一個空矩陣。空矩陣並不等於零,而是資料的另一種型式,可以用isempty()這個函數指令來檢驗。在輸入過程中,有時也必須利用這種技巧與程式使用者作順暢的交談。比如,輸入時常使用預設值,若使用者不輸入即直接按下enter,表示接受預設值,則程式可利用所得之空矩陣作檢驗,並設法餵給預設值。例如:

reply=input('要繼續嗎?y/n [y]','s');
if isempty(reply)
  reply='y'; %若為空矩陣,將其reply設為'y',以利程式進行
end


在input指令中,參數 's',表示輸入之內容當作字串,存於左邊的變數中。若沒有這個選項,則輸入之內容視為數值,此時若輸入文字則會有錯誤的信息。

7.2 MENU函數

menu這個指令為新版的功能。可改善input所無法顯示不同選擇的缺點。利用menu,可以將要選擇的項目逐一列出,並且先給一個標題或提示。經過選擇後,其值將轉存於左邊之變數中,其型式如下:

k=menu('標題', '選擇一','選擇二',…,'選擇N')


例如:

k=menu('請選擇你最喜歡的顏色:','紅色','黃色','藍色','白色','綠色')

















k = 5

這個指令適用於選擇性的問題,而且會自動產生一個小視窗,可以利用滑鼠選擇。選後得到的對應數值就存入k中(本例中選擇綠色,故k值為5),如此再作進一步的運算。

7.3 KEYBOARD & RETURN

keyboard是一個不用參數的指令。當程式執行到這個指令時,會將執行權開放到鍵盤,由鍵盤輸入。因此你可在鍵盤上下另外指令,查看資料、改變參數或檢查執行進度等等。等到在鍵盤上鍵入return指令時,執行權返回程式,在新更改的參數內容下,繼續往下執行。

7.4 DISP(X)

這個指令是最簡單列印資料的用法,其結果在螢幕中顯示。其中X可為矩陣,也可為文字模式,但不能兩者兼具。其顯示的結果不包含變數名稱,只有內容。例如:


>>disp('本期樂透彩中獎號碼為:')
ans=本期樂透彩中獎號碼為:
>>disp(rand(5,4)*10^3)
ans=
582.7917 225.9499 209.0694 567.8287
423.4963 579.8069 379.8184 794.2107
515.5118 760.3650 783.3286 59.1826
333.9515 529.8231 680.8458 602.8691
432.9066 640.5265 461.0951 50.2688


這樣的處理方式有點笨拙,文字顯示與矩陣顯示必須執行兩次,這是因為disp(X)中的X僅能處理一個參數。雖然如此,這個參數仍然可使用向量替換,故妥善的辦法是利用向量將一些數值經過轉換後,組成向量之型式以中括號起來。例如:
>disp(['今天的日期是 ' date]) 
今天的日期是 02-Apr-2006


若要文字與其他數值型變數共用,則必須利用num2str()的函數轉換;將數值轉成字串,然後用中括號括起來,當作一個向量參數處理。例如我的身高178公分,則:

>>height=178;
>>disp(['我現在的身高是 ' num2str(height) '公分。'])
我現在的身高是 178公分。


雖然經過這樣處理比較麻煩,但在程式寫作時,只要稍加用心也可以克服一些困難。這種處理方式同樣可以應用在input指令之文字串中,以提供更彈性的資訊給回答者。

7.6 列印函數FPRINTF

其格式如下:

count=fprintf(fid, format,A,…)


這是一個很通用的輸出指令,其輸出之格式依照指定的format,後面是所要輸出之參數名稱。其中fid代表輸出的裝置,fid=1代表標準裝置,或為螢幕,所以可以印在指令窗內。若要輸出至某一個檔案,則必須先將該檔案打開,即open指令打開一個檔,然後訂一個共通的代號fid。

列印格式則由format制定。format係由一系列包括C語言轉換碼之文字。其轉換控制包括符號、對齊、有效位數、欄寬及其他輸出之控制格式。此字串可包括特定字母,可代表非文字之控制標如新行及定位等。

這些轉換格式均以%符號開頭,並可能包括下列之項目:

  •  正負旗標
  •  欄寬及位數
  •  上下標
  •  轉換文字

上述的項目中,轉換文字必須存在,其餘視情況作選擇。












正負旗標位於輸出數字之前,可用+,-,0作為控制數字符號之用,若使用負號(-),則表示採用左邊對齊,正號(+)代表要印出該數字之正負號,若用零(0)置於前,表示其前面滿0。其例如%-10.5d、%+10.5d、%010.5d等等。

欄位寬度及有效位數則由數值表示,小數點前之數字代表總位數,如%6f,表示要列印六位數;若小數點後面有數字,表示該數之小數位數,如%6.3f表示有三位數小數,但總欄寬(包括小數點)為六位數。

格式標則在標明印出數值之性質,常用的如d代表十進位,e代表指數表示法,相關資訊如表5.2:

表7.2 格式標代碼表














格式標代碼 說明
%c單一文字
%d十進位(有正負符號)
%e 指數格式(如小寫1.414e+01)
%E 指數格式(如大寫1.414E+01)
%f 固定位數
%g 較%f或%e更為簡潔的型式,使用小寫e
%G 同%g但使用大寫E
%i 十進位(有正負符號)
%o 八進位(無正負符號)
%s 一串文字
%u 十進位(無正負符號)
%x 十六進位(使用小寫文字a-f)
%X 十六進位(使用大寫文字A-F)


在format格式中,也常使用控制標以印出控制的動作,這些包括換行、換頁等等,其項目如表5.3:

表7.3 控制標代碼表









控制標說明
\b倒退一格
\f進一頁
\n跳一新行
\r回行
\t水平跳格
\\反斜線符號
\''單括號
%%百分比%


下面為產生一串正弦值,然後將其存於一個名叫sinx.txt之文字檔案內。在新開一個檔時,通常要先用fopen開啟,等利用fprintf存完後,再用fclose令其關閉。才算存檔完成。

x=0:pi/10:pi; % %將pi分成十等分,置於x中
y=sin(x); %計算sin(x)值,置於y
% 打開一個名叫sinx.txt的文字檔,將其代碼設為fid,'w'代表要寫入
fid=fopen('sinx.txt','w');
%列印資料x,y於檔案中,格式x為%6.4f;y為%10.8f
fprintf(fid,'%6.4f %10.8f\n',x,y);
fclose(fid); %關閉檔案

type sinx.txt %將存檔內容印出

0.0000 0.31415927
0.6283 0.94247780
1.2566 1.57079633
1.8850 2.19911486
2.5133 2.82743339
3.1416 0.00000000
0.3090 0.58778525
0.8090 0.95105652
1.0000 0.95105652
0.8090 0.58778525
0.3090 0.00000000


在format格式中,亦可包括其他文字,使印出的報告容易閱讀。比如說:

>>height=178;
>>fprintf(1, '我的身高有 %4.0f cm,也就相當於 %4.2f ft.', height, height/30.48);
我的身高有 178 cm,也就相當於 5.84 ft.


這樣的文字讀起來也比較人性化一些。在此,fprintf括號內第一個參數1表示印在螢幕中,而中間即為格式之內容。其中除文字本身外,遇到%時表示是輸出變數的格式,而每遇到一個%的型式就需要有一個變數對應。故第一個%4.0f 對應變數height,第二個%4.2f對應height/30.48。

由於格式format 中,已經用掉撇號'、反斜線\及百分比號%,因此若要在文字中出現這些文字時,必須寫上兩次。例如:


>>fprintf(1,'My mother''s book is not 100%% complete yet.\n')
My mother's book is not 100% complete yet.


最後的\n表示跳到下一行列印。
相對於輸出之fprint指令,讀入資料可以利用fscan指令,其應用請待下回分解。

7.7 fscanf讀取資料

正式之讀取資料函數如下之格式:

A=fscanf(fid,format)

[A, count]=fscanf(fid,format,size)

這個指令可以自fid所標示的檔案中將其資料依format的格式取出,並置於矩陣A中。fid的定義與前述之fprintf指令相同。另一種型式則增加count與size兩參數。count表示已完成的資料筆數。而size則是決定讀入之資料量,可用[m,n]表示,表示讀入之資料可填滿m x n的矩陣。其中n 可用inf取代,代表EOF(End of File)。

此處格式format%符號大體上與fprintf相同,其中%e,%f,%g均代表具有浮點之數據。以%12hd為例,12為數值之總位數,d代表使用整數,其前面所置的文字表示整數的型式,h代表短整數,l代表長整數,而lg則代表雙精度浮點值。這些文字代碼可參考fprintf指令,其功能相似。


fid=fopen('sinx.txt');
A=fscanf(fid,'%g %g',[2 inf]); %最後項為size,表示讀入二列資料,直到檔案底
fclose(fid);

A=A'
A =
0 0.3142
0.6283 0.9425
1.2566 1.5708
1.8850 2.1991
2.5133 2.8274
3.1416 0
0.3090 0.5878
0.8090 0.9511
1.0000 0.9511
0.8090 0.5878
0.3090 0

7.8 存取excell檔案資料

許多資料多以微軟的速算表excell儲存,故如何自速算表中取得資料相當重要。為此,MATLAB有如下之指令可以自速算表中擷取資料:


[N, T, rawdata] = xlsread('filename', sheet, 'range')


其中filename為檔案名稱,包括路徑及副檔名.xls;sheet為速算表之序號或名稱,若為名稱則必須使用字串型式,如'sheet1'是。range則為所要指定的範圍,如'A1:C3',必須為範圍的型式,大小寫均可。若不指定範圍,即表示輸入整個速算表。若sheet=-1,則有特別功能,亦即可以直接在速算表指定其範圍即可,但此時必須將速算表同時開啟才行。在輸出參數方面,N為數值陣列,T為字串陣列,而rawdata則為細胞陣列,包括所有資料。

下面舉一個例子。設一速算表之內容如下:






















>> [N,A,rawdata]=xlsread('demo7.xls')
N =
12 NaN 60
14 NaN 105
16 NaN 56
18 NaN 48
20 NaN 50
A =
'頁數' '姓名' '年齡' '說明'
'' '馮丁樹' '' '本書的作者'
'' '張三豐' '' '太極拳祖爺'
'' '曹雪芹' '' '紅樓夢作者'
'' '高山青' '' '一首常唱的歌'
'' '文天祥' '' '作了正氣歌'
rawdata =
'頁數' '姓名' '年齡' '說明'
[ 12] '馮丁樹' [ 60] '本書的作者'
[ 14] '張三豐' [ 105] '太極拳祖爺'
[ 16] '曹雪芹' [ 56] '紅樓夢作者'
[ 18] '高山青' [ 48] '一首常唱的歌'
[ 20] '文天祥' [ 50] '作了正氣歌'


其中,N為數值矩陣,A為字串矩陣,而rawdata則為全部之資料,亦即細胞矩陣。


輸出至檔案之格式與參數也類似讀入之型式,其指令格式如下:

[status, message] = xlswrite('filename', M, sheet, 'range')

除上述之定義外,其中M為待輸出之矩陣,輸出檔案名稱可包含路徑。左邊有兩項資訊,若輸出成功,則status=1否則為0,message則會傳出錯誤的信息。其中有兩項,一為錯誤資訊的內容,另一為identifier,或為識別項。通常在寫入excell檔時,若原檔已有內容,則在進行存檔後會有如下之警告信息:

Warning: Added specified worksheet.

若不想這類信息一再出現,可以打入如下指令:

>>warning off MATLAB:xlswrite:AddSheet

如此,警告信息不會再出現,而後面MATLAB:xlswrite:AddSheet即為識別項idengifier。


例如:


>> books = {'書名', '價格(元)'; '模貝設計' 298; '模貝設計' 298; '基本電學' 580; '哈利波特' 660}
books =
'書名' '價格(元)'
'模貝設計' [ 298]
'模貝設計' [ 298]
'基本電學' [ 580]
'哈利波特' [ 660]

>> status=xlswrite('demo8.xls',books,'demo')
status = 1 %表示存檔成功


在excell之檔案內容如下:

10/25/2006

第六章 自訂函數

6.1敘述檔與函數檔


MATLAB依執行時之實質型式分為兩種,其一為前面所述之敘述檔( script files),另一為函數檔(function files);兩者均可在指令窗中呼叫。但敘述檔中,使用者可將檔案內所含的指令一一地在指令窗下執行,其執行過程中,所發生的變數均為整體性變數。因此執行敘述檔時,只能呼叫其檔案名稱,不能在名稱後面加上任何輸入參數。若要事先設定一些參數初值,辦法之一是執行之初立即設定,例如在工作空間內更動其中之變數值,然後再執行該指令。其二是在檔案中利用指令讀取特定檔案之資料或由鍵盤輸入資料。由鍵盤輸入指令可以累積成為指令之組合,並給予指令群之名稱,此名稱因而亦可作為指令在指令窗中直接執行。例如,有一個敘述檔之名稱為 pick_a _number.m,則在 MATLAB 指令窗中直接鍵入"pick_a _number"即可執行該檔案之內容。


% pick_a_number.m, a script file for demonstration.
% Pick up a number that equals to the randomly generated one
%
message='Please pick a number from 0-9 ==> ';
n=1;
while 1
A=input(message);
if isempty(A), A=0; end;
pick=fix(rand*10);
if A==pick,
disp('Congradulations, it''s done!')
break
else
disp(['The random number is ', num2str(pick),...
'. You have failed ',num2str(n),' times!'])
n=n+1;
end
end


執行的結果如下:

>> pick_a_number
Please pick a number from 0-9 ==> 2
The random number is 6. You have failed 1 times!
Please pick a number from 0-9 ==> 3
The random number is 2. You have failed 2 times!
Please pick a number from 0-9 ==> 4
Congradulations, it's done!


上述的程式為一個敘述檔,其內容為讀者自訂一個數字,由程式隨機產生一個數字,直到兩個數字相同為止。此程式使用無窮迴圈while 1,…end,故除非答對,才利用break中斷跳出,否則會繼續進行。有關迴圈指令,將在另節詳加討論。

在敘述檔中之變數是整體性(global) 的,也是共用的 ,故任何先前定義的變數,均可為敘述檔執行時所用,或改變其內容,所以它不需考慮輸入及輸出之問題。敘述檔用於必須重覆多次使用一群指令時最為方便。

函數檔案的內容與敘述檔大略相同,但其檔案可以作為MATLAB語言之一部份,可有輸入及輸出,且其變數的空間是自主的,不與公共空間共用。在型式上函數檔案之開頭一行需有一個函數名稱之宣告,然後用小括符包括其需要之輸入參數。函數名稱前可以有等號,在此等號之左方為其輸出之參數。等號之右方為函數名稱,其後有輸入參數。這些參數變數,無論是輸入或輸出,均可能代表一個數值或一個矩陣。在其他程式語言中,矩陣之變數均需經過宣告,在MATLAB的世界裡,可以省去此項運作,但若屬於大矩陣,事先安排矩陣空間,對提高運算效率有相當的助益。輸出入之參數中,均會預設為矩陣的型式。在函數檔案中,其所屬之變數均屬區域性的,除非特別另宣告為整體性參數。相關的函數特性等後段再作論述。

函數檔案之類型純為文字檔,故可在任何文書處理軟體中編輯或修改。不過,MATLAB也有一個程式編輯器,可以在此編輯器中偵錯(Debug)或編輯,這個編輯器比其他文書處理具有特殊的功能,新版中並有行號。若不使用編輯器,也可以在指令窗中亦直接打入type之指令,以觀察該檔案之內容。

6.2 自訂函數之型式

在MATLAB中,可以輕而易舉地定義自己需要之函數,這點大大增強MATLAB之應用能力。一個自定的函數為一個M-檔案,其儲存名稱須與函數名稱相同,格式如下:

function [輸出變數]=Name_of_function(輸入參數)

函數的名稱需以英文字母開頭,中間可為數字或底線,但其間不能有+,-,*,/等字眼。函數名稱長度依一般檔案名稱之規定,但最多僅認定63個字,不過長度如此長也不盡實際。且其名稱不要與現有MATLAB內定的相同,通常可用isvarname這個函數先檢查。輸出及輸入參數可為多個,亦可為矩陣,但輸入參數必須用一般括號括起來,而輸出參數若僅有一項輸出,則不必用任何括號,若有兩項以上,則需用中括號括起來。下面為一個簡單的例子:


function y = freebody(time)
% calculation of height with respect to time
% time: the elapsed time in seconds
% y:the depth of falling, cm
% Example: yy=freebody(1:20)
y =1/2*980*time.*time;%y=(1/2)gt^2




這是一個計算自由落體的例子。freebody為函數之名稱,time為其輸入變數,可為矩陣或向量。這個函數以time=[0:4:20],即time=[0 4 8 12 16 20]代入,執行如下:


>> yy=freebody(0:4:20)

yy =

0 7840 31360 70560 125440 196000



特別注意的是四秒時,自由落體將降至78.4m;若為4.5秒時,約達100m。911時紐約世貿大樓幾乎是以這個自由落體的速度塌下來,沒有任何阻檔,所以一直是個謎。在函數內容中,其說明行之前以%開頭,其後面的資料程式即認為是說明,不加執行。依照慣例,函數除第一行外,其最開始不包括空行的幾行為說明行,但每行前面需加%符號。這個說明行主要在敘述該函數之功能及使用方法。當你使用help freebody查詢這個自設的函數指令時,會出現在其說明之內容。例如:


>> help freebody
time: the elapsed time in seconds
y:the depth of falling, cm
Example: yy=freebody(1:20)




程式本體則包括其他函數呼叫、程式結構流程及交談式輸出入、計算、設定指令、註解及空行等。註解可佔全行,亦可在敘述指令之前後。但其前面均需加上%。MATLAB發現有%存在時,其後面之敘述均視為註解。若說明行不需佔有全行,則亦可在指令後加%後,隨時作進一步的說明,其情況如下:


   y =1/2*980*time.*time; %y=(1/2)gt^2



若有多行的註解行,但不想每行使用%為開頭時,可以在註解之前後行加上%{及%},這兩個符號必須單獨佔行。此外,MATLAB亦不處理空行,所以為使程式容易閱讀,可以加空行。但在程式開頭的第一次空行會中斷help指令所顯示的內容。

函數中所用之變數若不在函數之參數行列出,則應屬區域性變數。若要成為整體性變數則需要特別宣告(如global x y)。下面為兩連桿相接之位置的計算,在此自訂函數two_vectors中,其輸入項為r1,r212,輸出為結果點B之位置(x, y)。程式中,另有自創變數d2r, th1, th2等,這些均為區域性變數,故在運算結束跳離函數之後,這些變數將與函數外之參數無關。






function [x,y] = two_vectors(r1,r2,theta1,theta2)
% Find the resultant position of two vectors
% r1,r2: lengths of the two vectors, cm
% theta1, theat2: horizontal angles of r1 & r2, deg.
% x,y:coordinates of the tip position.cm
% Example: [xx,yy]=two_vectors(5,8,30,60)
d2r=pi/180;
th1=theta1*d2r;th2=theta2*d2r;
x =r1.*cos(th1)+r2.*cos(th2);
y =r1.*sin(th1)+r2.*sin(th2);



上項例子設兩桿之長度為5cm及8cm,對應水平角為20及50度,則經執行指令之結果如下:


>> [xx,yy]=two_vectors(5,8,20,50)
xx = 9.8408
yy = 7.8385


6.3 函數之編輯

前面所設之函數名稱如freebody、two_vectors等在存檔時,必須直接以函數名稱為其檔案名。以後呼叫時才能找到這個檔案來執行。存檔時,必須以這兩檔的名稱加上”.m”的字尾,即如freebody.m、two_vectors.m。撰寫這些函數檔可用MATLAB的編輯器(editor),只要從檔案指令下開檔即可找到。或直接在指令窗下這樣的指令:

>>edit two_vectors.m


編輯器可同時開啟許多檔,但僅能執行其中一個檔。編輯窗具有行號,可以進行偵錯及修改,相當方便。作為註解的部份則會以不同的顏色顯示,只有真正程式內容部份用黑色(預設值),但其行號仍然由第一行算起。在偵錯期間,此行號成為指向錯誤點的定位方法。編修完後的檔案即可存成.m檔。

進行偵錯時,可以採用單步法或全程執行模式,但可設中斷點。只要用滑鼠到行號右側一點,該會顯示一紅點。將來程式執行時,抵達中斷點會暫時停止。此時可以利用鍵盤下各種指令以確定執行過程中各項變數的變化。中斷點可以設好幾點,執行時會依序暫停。要解除中斷點,只要重複用滑鼠再選一次即可,或指向該行然後在指令行上按打叉的指令即可。單步法選下圖第三項,每按一次僅執行一行,如此可以偵錯全部。第四項與第三項相同,但可進行呼叫的函數內部逐步執行。第六項為往前執行,直接終了,或遇到紅點的行為止。其餘讀者可以自行試驗其功能。

MATLAB的新功能中,有一項是可以檢驗所選寫程式的執行效率。這個功能置於編輯器的工具項中的一個選項,稱為M-LINT檢驗程式碼(Check code with M-Lint)。這個選項可以在程式完成偵錯後,再進一步檢查程式的執行效率。由此可以得到一個建議的報告(M-Lint Code Checker Report),其內容有相關修正事項,並如何使程式執行時更有效率,且節省處理時間等等。對於大程式而言,這項功能相當有用。

6.4次函數

一個函數檔案可以包含一個以上的函數,這些多餘的函數又稱為次函數或副函數(Subfunction)。這些副函數僅能為主函數或其他同檔案之次函數相互呼叫。每個次函數均有其函數之定義行,而且一個接一個置於主函數後面,其順序不拘。下面為一個畫圓的例子,輸入圓的半徑,然後以亂數函數決定其圓心位置,再利用plot指令繪出圓。程式先從drawcircles(r)進入,在其for…end之迴圈中呼叫次函數circ(rr)決定圓心位置及計算圓的座標,而在其中又呼叫randxy之次函數,計算圓心位置,此時即為次函數呼叫次函數之例子。在主函數中利用plot指令則只要將對應之X、Y座標值輸入就可,其結果如圖。


function drawcircles(r) % 主函數
% Draw circles with radii of r.
global MAXR;
n=length(r);
MAXR=max(r);
hold on;
for i=1:n,
[X,Y]=circ(r(i));
plot(X,Y);
end
hold off;
axis equal;
end

function [x,y]=circ(rr)%次函數
% Calculate the points of a circle.
[xx,yy]=randxy;
theta=linspace(0,2*pi,60);
x=xx+rr*cos(theta);
y=yy+rr*sin(theta);
end

function [xx,yy]=randxy%次函數
global MAXR
% locate the position of the center.
xx=rand*MAXR;
yy=rand*MAXR;
end

執行此程式:

drawcircles(2:10)






















由於使用亂數決定各圓之圓心,故每次執行上項程式均會有不同的結果。只是其半徑均會一樣,本例之半徑為[4 3 1 5]。程式中,除兩個次函數circ 與randxy外,尚有呼叫max及length兩函數,這些都是MATLAB的內在函數。

次函數內之變數必須透過輸出入參數傳遞,否則僅能在自家之函數範圍內使用,包括主函數及其他次函數均無法讀取該次函數內之變數。除非將要共用之變數宣告為全域性(Global)。例如,函數中有用到全域變數MAXR,必須在使用到的函數(如randxy)中要宣告,才能共同使用。

呼叫次函數時,MATLAB首先檢查是否為次函數,然後檢查具有相同名稱之自有函數,最後才是內建函數。故若有相同函數名稱仍以次函數為優先。

6.5匿名函數

匿名函數是一種不必具有函數檔的型式即可進行呼叫的函數。因此建立匿名函數可以在敘述檔或函數檔中加入,甚至在指令窗內加入即可。其基本型式如:

fhandle=@(arglist) expression

這個型式右邊以@開始,利用這個符號產生函數之握把(handle),以此作為呼叫此函數之起始,形成匿名函數之一部份。括號內之arglist為自變數名單,是由外進入函數之變數。而經空一格後之expression則為函數之主體。此匿名函數之握把可以設定給左邊之變數fhandle,利用此變數呼叫此匿名函數。下面為例子:


>>circ=@(r) r.*r*pi;
>>circ(2:2:10)
ans =
12.5664 50.2655 113.0973 201.0619 314.1593


這是一個計算圓面積的例子,只要輸入半徑,即可得到圓的面積。這時之自變數為r,而函數握把為circ,也是匿名函數之名稱。利用此握把,即可下指令執行該函數,因此直接將其作為函數之名稱餵入參數項即可執行。例如計算一個圓環所佔的面積,設內環半徑為2、外環為4,則以LoopArea可得:


>>LoopArea=circ(4)-circ(2)
>>LoopArea =
37.6991


若為兩個自變數,則其變數可以酌增,如:


>>Cllipse=@(x, y) 3*x.*x + 5*y.*y;
>>Cllipse([3 4 6],[5 7 8])
ans =
152 293 428


利用匿名函數,以小型函數為主,故使用上相當方便。這種函數亦可分組在細胞陣列之中,隨時呼叫所需進行之計算,例如:


>>Afun={@(x) x.*x, @(y) y.*y+100 , @(x,y) x.*x + y.*y + 100}
Afun =
Columns 1 through 2
[@(x) x.*x] [@(y) y.*y+100]
Column 3
[@(x,y) x.*x + y.*y + 100]
>>Afun{1}(5) + Afun{2}(10)
ans = 225

>>Afun{3}(5,10)
ans = 225


由於所設的匿名函數值第三項等於第一及第二項之和,故上述所求得之結果應會相同。

匿名函數定義時,其內也可以使用常數。但這些常數必須事先定義,才能完成匿名函數的型式。此匿名函數之常數將維持原先設定之值,不再改變。故若在運算過程中間即使更改該常數值,也不更動匿名函數之原常數值,除非重新定義該匿名函數。以下為例:


>>a=2;b=3; f=@(x) (a*x+b)
>>f = @(x) (a*x+b)
>>f(2)
ans = 7

>>a=5;b=8;f(2)
ans = 7

顯然這個匿名函數並未因常數a與b之改變而改變。必須必須在常數改變之後重新定義才能生效:

>>a=5;b=8; f=@(x) (a*x+b)
f = @(x) (a*x+b)

>>f(2)
ans = 18


在MATLAB有繪圖指令fplot可以繪製一般函數,其輸入只要函數及其前後之範圍:


>>a=3;b=4;g=@(x) (a*x.^2 +b*x)
g = @(x) (a*x.^2 +b*x)
>>fplot(g, [-10 10])





















也可以直接將匿名函數整個內容置入於其第一參數之中:

>>fplot(@(x) (a*x.^2 + b*x), [-20 20])


6.6巢狀函數

在函數中又存在另一個函數,這是MATLAB函數中基本的型態,不過大部份是另一個函數附於主函數之下,因此可以隨時呼叫。但此處所謂巢狀函數(Nested function)是處於函數內文之中,或以自己為函數重新呼叫。

撰寫一個巢狀函數與一般函數相同,亦即在函數中鉗置一個格式相同的函數:


function mainfun=A(p1, p2)%主函數
E(p1); B(p1); D(p2);A(p1,p2)
function infun1=B(p3) %第一層巢函數
C(p3);E(p3);B(p3);D(p3);A(p3,p3);

function infun1=C(p4) %第二層巢函數
     C(p4);E(p4);B(p4);D(p4);A(p4);
    end
end
function infun1=D(p3) %第一層巢函數
  F(p3);E(p3);B(p3);D(p3);

function infun1=F(p4) %第二層巢函數
     F(p4);E(p4);B(p4);D(p4);A(p4,p4);
    end
end
---
end

function main2=E(p5) %次函數
E(p5); A(p5,p5);
end


一般函數並不需要在終了加一個end,但若要做成巢狀結構,則必須以end結束。中間之函數可以再包含另一個巢狀函數,或者最外圍亦可包含多個巢狀函數。

在呼叫層次則有限制。第一層主函數可以呼叫次函數E及其下之第二層巢B與D函數,但不能呼叫其下之第三層以下如C函數。低層函數則可以呼叫其本身及同層以上之高層函數(如D可呼叫D、A、B、E等函數。

函數中之變數影響範圍則依其巢函數之位置而定。通常主函數與次函數各享有自己的工作空間,故除非利用輸入參數傳介,或宣告全域性變數,否則這兩個工作空間是相互獨立的。巢函數則與次函數不同,因為它與其直屬函數共享有工作空間。因此,巢函數本身雖有自己的獨立的空間,但由直屬函數內之變數則是共用的,因此這些變數均可能在兩個函數裡被更改。

巢函數置於某一函數之中,在執行時除非經由呼叫程序,否則不會自動執行。因為它也具有閉鎖性質,其變數也需經由輸出入參數傳輸。下面為一個計算多項式值的函數nest_fun。其輸入x 值可以為列矩陣,而係數a有三項,分別為A、B、C ,a可以有多項輸入。如此計算可以得到不同的組合。本函數利用細胞陣列rr_arry作為儲存運算結果,讀者可試著瞭解前節所介紹的功能。要顯示其內容必須使用celldisp指令。


function [rr_array]=nest_fun(x,a)
%function to find sets of polynormials.
% a: set of constants, [A B C]
% x: variables in array
% Example: rr=nest_fun(2:10,[1 2 4;2 4 8])
n=size(a);
for i=1:n
A=a(i,1);B=a(i,2);C=a(i,3);
rr_array{1,i}=['A=',num2str(A),', B=',...
num2str(B),', C=',num2str(C)];
rr_array{2,i}=polyx(x);
 end
function r=polyx(xx)
  r=A.*x.^2 + B.*x +C;
end
end


執行例:


>> rr=nest_fun(2:10,[1 2 4;2 4 8])

rr =
'A=1, B=2, C=4' 'A=2, B=4, C=8'
[1x9 double] [1x9 double]

>> celldisp(rr)

rr{1,1} =A=1, B=2, C=4
rr{2,1} = 12 19 28 39 52 67 84 103 124
rr{1,2} =A=2, B=4, C=8
rr{2,2} = 24 38 56 78 104 134 168 206 248

6.7函數輸出入參數

在函數應用中,變數之轉換至為重要。一般函數內定義之變數均屬區域性變數(Local variables),只能提供該函數內部使用,當函數執行完畢後,即被清除。為維持參數之一貫性,當呼叫某函數,也必須同時傳遞重要之變數作為該函數之輸入值,以為繼續計算之根據。執行函數完畢時,有些結果仍寄望能傳回呼叫主程式,此時必須有輸出參數。

6.7.1全域性變數


有些變數對許多函數而言有時同等重要,此時就可宣告為全域性變數(Global variables),並在需要此變數的函數中同時宣告。此變數即可共同使用,不會因函數之轉換而遭到清除。在前面所討論的敘述檔中,實際上其所使用的變數均為全域性的變數。因此每次執行後,其相關之變數均有可能被更改。

全域性變數之宣告如下例:

global GRAVITY X Y Z


上式即宣告GRAVITY、X 、Y、Z等四個變數為全域性,習慣上均用大寫,以與其他區域型變數區別。不過,這並不是必然的作法,小寫也一樣可用。

6.7.2 eval函數


一般下指令時,通常以其指令名稱及需要之參數在指令窗下或在檔案中直接下,但有些時候我們希望指令也是一種變數的型式,使指令內容也能依檔案之內容而變化。換言之,先將指令化為文字串,但然依文字串之內容轉為指令去執行。MATLAB提供一個eval函數具有轉換的功能,其格式如下:

eval(指令相關敘述)

例如:以下面程式求某一連續正整數累進之總和


% evalsum.m for demonstration
% Find the sum of cumulative sum(1:n)
tsum=zeros(1,10);xsum=0;
for i=1:10,
xsum=xsum+eval(['sum(1:',int2str(i),')']);
tsum(i)=xsum;
end
1:10
tsum


執行後:


>> evalsum
ans = 1 2 3 4 5 6 7 8 9 10
tsum = 1 4 10 20 35 56 84 120 165 220


式中之eval函數是將sum(1:1)、sum(1:2),…sum(1:10)逐一計算後累加。由於其中sum之參數逐一變化,故必須用數值改變為字串後(int2str),再用eval函數轉為指令執行之。

6.7.3 輸出入參數個數nargin與nargout


呼叫函數時,大部份需要輸入參數,經過處理過後,也需要輸出參數以獲得希望的結果。為使程式能掌握呼叫的參數,通常程式中有設定輸入與輸出參數之個數。其變數分別為nargin與narout,使程序進行前能依實際需要,決定所要進行的步驟。所以上述兩參數名稱已為程式內部參數,程式中可依需要使用,但在設定變數時,需不得重覆。

下面的例子為計算空心圓柱體體積之程式,主要針對輸入項少於原設值時,採用預設值的方式。


function [volume]=pillar(Do,Di,height)
% Find the volume of a hollow pillar
% ro,ri:outside & inside diameters
if nargin<1,
height=1;
di=0;
do=1;
end
volume=abs(Do.^2-Di.^2).*height*pi/4;


此程式利用nargin變數,確定所缺少的輸入項並補上預設值。所以即使不使用輸入參數,其結果將以直徑及高度均為一單位之圓柱體計算其體積。如:


>>pillar
ans =
0.7854
>>pillar(5,2,[1:5])
ans =
16.4934 32.9867 49.4801 65.9734 82.4668


即使利用pillar函數,亦可計算陣列輸入時之對應體積。可以增加程式之使用彈性。對於nargout參數之應用亦與nargin相同,當下指令時,有時未提出輸出的項目,此時其nargout=0。此時亦可利用這種情況決定是否進行繪圖或作其他動作。

6.8私有目錄函數

在MATLAB的目錄中,有一個private子目錄,任何函數置於此目錄中,可以獲得優先執行。由母目錄可以事先看到其中的函數,因此主函數也可以利用此目錄之區隔獲得執行之優先權。為此若有兩相同名稱之m-檔案,一若置於private子目錄而另一置於其他目錄。則在private子目錄則在下之同名稱檔案會優先被執行,置於其他目錄檔案者雖然名稱相同,亦不會被執行。

private子目錄通常無法用檔案總管觀察,但可以用下列指令得到訊息:

>>help private/myprivfun

6.9函數為輸入參數之呼叫法

有些函數之輸入參數必須呼叫其他函數名稱,以該函數進行運算。例如fzero、fminsearch、 fminbnd等等,其中均需輸入或呼叫適當函數名稱,使運算正常進行。大體上傳輸函數之方法有四:

一、 以字串之型式對照M函數檔之名稱。例如:


function y = freebody(time)
% calculation of height with respect to time
y =1/2*980*time.*time;

>>[x,f]=fzero('freebody',[0 1])
x = 0
f = 0


二、 以現存函數名稱之握把呼叫,例如:


>>[x,f]=fzero(@freebody,[0 1])
x = 0
f = 0

>>f1=@freebody;[x,f]=fzero(f1,[0 1])
x = 0
f = 0


三、 利用inline指令轉換,例如:


  >> g=inline('freebody(time)')
g = Inline function:
g(time) = freebody(time)

>> [x, value]=fzero(g,[0 1])
x = 0
value = 0


四、 直接使用運算式,例如:

>>f1='x.^2-5';[x, value]=fzero(f1,[0 5])

x = 2.2361
value = 8.8818e-016

或,

>>[x, value]=fzero('x.^2-5',[0 5])
x = 2.2361
value = 8.8818e-016


上述四種呼叫方式原則上均可適用。但在MATLAB之應用上,仍以握把呼叫較有效率,其執行速度也較快。雖然如此,由於撰寫程式本乎於習慣,故何者為佳,端視當時情況而定。一般之匿函數屬於握把呼叫型,可以使用在任何主函數內的位置。但這種函數以簡短為主,無法應用在大程式上。以下面之三個匿名函數為例,f1、f2各自定義內容,而g則由f2配合f1,輾轉以x為自變數,則不必利用傳統的函數型式,可以輕易加以串聯,並得到g([1 2 3])之結果:


>>f1=@(x) sqrt(x);
>>f2=@(x) 5*exp(-x);
>>g=@(x) f2(f1(x));
>>g([1:3])

ans =
1.8394 1.2156 0.8846

10/24/2006

函數之解說與示範

這是一組程式偵測某一點是否在圓圈內或在方塊內,並作判斷的示範例。

10/22/2006

第五章 資料結構

5.1資料型式


在MATLAB中,其基本資料型式約有十五種,每種均可以陣列與矩陣表示。有些前面已經使用過,有些則屬使用者定義的範圍。一般的資料型式除整數,有單精度、雙精度、邏輯值、字母、細胞陣列、結構、函數握把值等等,視需要而定。數值的表示法則有整數、浮點、複數等等。對於任意數值屬於何種型式則可利用class函數進行檢驗,在指令窗中亦可利用whos函數顯示其大小與型式。

5.2多維陣列

5.2Multi-Dimensional Arrays


一般的矩陣多探討兩維。三維陣列雖僅增加一個維數,但型式就複雜許多。二維陣列以行、列表示,故在標示上,僅有兩個下標,如A(3,4)即為列3行4的表示法。若增加一維,則必須再增一個下標位置。而其所代表的意義變成另一層,而不是行列的觀念,增加之層次通常以頁數(pages)表示之。故A(3,4,1)變為第一頁,而A(3,4,2)為第二頁。其前兩標仍維持二維之行列觀念,故A(3,4,1)表示第一頁第三列第四行之元素,而A(3,4,2)表示第二頁第三列第四行之元素。

若要指出第一頁或第二頁所有的元素,則可以使用冒號代表。如A(:,:,1)、A(:,:,2)是。

>>A(:,:,1)=[1 2 3;4 5 6;7 8 9]; %第一頁輸入,結果不顯示
>>A(:,:,2)=magic(3) %第二頁輸入,並顯示整個矩陣之結果
>>A(:,:,1) =
1 2 3
4 5 6
7 8 9

>>A(:,:,2) =
8 1 6
3 5 7
4 9 2

輸入第一頁時其頁數可以省略,處理其餘頁數時則必須標明。通常要觀察一個多維矩陣之維數,可用ndims或size函數,只是兩個結果之表示法略有不同,例如:

>>ndims(A)
ans = 3

>>size(A)
ans = 3 3 2

多維矩陣之形成,亦可利用另一個函數cat(dim, C, D,…)進行增減,cat又稱為串接指令,可以將許多小矩陣依dim維度之安排進行組合。dim參數即表示所要做的矩陣維數。當dim=1時,表示產生一個[C;D]疊加的行向矩陣;dim=2時,則產生一個[C D]疊加的列向矩陣;而dim=3時,則產生一個分別由C與 D各佔一頁的三維矩陣。例如:

>>c=[1 2;3 4], d=[5 6;7 8]
c =
1 2
3 4
d =
5 6
7 8

>>cat(1,c,d)
ans =
1 2
3 4
5 6
7 8

>>cat(2,c,d)
ans =
1 2 5 6
3 4 7 8


>>cat(3,c,d)
ans(:,:,1) =
1 2
3 4

ans(:,:,2) =
5 6
7 8


在組成多頁之多維系統時,基本上其大小應一致,否則無法搭配。以下面之列為例,在第三維中,其二維之大小必須相同:

>>A=cat(3,[9 2;6 5],[7 1;2 8])

A(:,:,1) =
9 2
6 5

>>A(:,:,2) =
7 1
2 8

>>A=cat(3,[9 2 4;6 5 3],[7 1 2;2 8 4])
A(:,:,1) =
9 2 4
6 5 3

A(:,:,2) =
7 1 2
2 8 4

>>B=cat(3,[3 5;0 1; 1 4],[5 6;2 1;2 8])
B(:,:,1) =
3 5
0 1
1 4

>>B(:,:,2) =
5 6
2 1
2 8


若三維仍然不足,必須增加為四維時,則必須增加一個表示第四維之下標。在置放時,其第三維除非特別說別,均設為1,而其餘矩陣則依第四維之順序置入,亦即所輸入之資料是置於第三維之第一頁中:

>>C=cat(4,[4 2;6 5],[5 8;8 2])
C(:,:,1,1) =
4 2
6 5

C(:,:,1,2) =
5 8
8 2

>>C=cat(4,[1 2 5;4 5 7],[7 8 0;3 2 9])
C(:,:,1,1) =
1 2 5
4 5 7

>>C(:,:,1,2) =
7 8 0
3 2 9

若M與N為兩三維矩陣,今擬置於另一四維矩陣中時,其結果如下:

>>D=cat(4,A,B,cat(3,[1 2;3 4],[4 3;2 1]))
>> A=cat(3,[1 1;2 2],[3 3 ; 4 4])
A(:,:,1) =
1 1
2 2
A(:,:,2) =
3 3
4 4
>> B=cat(3,[11 22;33 44],[44 44;55 55])
B(:,:,1) =
11 22
33 44
B(:,:,2) =
44 44
55 55
>> D=cat(4,A,B,cat(3,[7 7; 8 8],[ 9 9; 2 2]))
D(:,:,1,1) =
1 1
2 2
D(:,:,2,1) =
3 3
4 4
D(:,:,1,2) =
11 22
33 44
D(:,:,2,2) =
44 44
55 55
D(:,:,1,3) =
7 7
8 8
D(:,:,2,3) =
9 9
2 2


其中,D(:,:,:,1)是置放M矩陣;D(:,:,:,2)是置放N矩陣;而D(:,:,:,3)則置放後面所串接之矩陣,其維數可以檢測如下:

>>ndims(C) = 4
ans = 4

要取出其中之陣列應用時,必須在整數的指標指到其對應之陣列,例如:第四層中之第二頁、第三層中之第二頁之第一列第二行為:

>>D(1,2,2,2)
ans = 6

>>D(1,2,2,3)
ans = 3

>>D(1,2,1,2)
ans = 5

為方便起見,超過二維之多維部份亦可用指標向量指向其頁數矩陣:

>>k=D(:,:,1,[1 3])
k(:,:,1,1) =
9 2
6 5

k(:,:,1,3) =
1 2
3 4

5.3細胞陣列

5.3Cell Array


細胞陣列與其他矩陣之情況略有不同,它在同一個陣列名稱下,儲存不同類型的資料陣列,諸如字串、文字、數值、複數、正整數陣列等,由於陣列之維數也容許不同,其應用的範圍更具彈性。有關細胞陣列之指令如下表:












指令型式 說明
A=cell(n) 產生一個n x n的空細胞陣列A
A=cell(m,n) 產生一個m x n的空細胞陣列A
celldisp(A) 顯示細胞陣列A之內容
cellplot(A) 繪出細胞陣列之佈置圖
[x,y,…]=deal(A,B,..) 進行串列配對,即x=A,y=B
[x,y,…]=deal(A) 進行串列配對,即x=A,y=A
iscell(A) 檢查A是否為細胞陣列,若是則傳回1,否則傳回0
A=num2cell(B) 將一數值矩陣B轉為細胞陣列A
cell2struct 將細胞陣列轉換為結構陣列
cellfun 應用細胞函數於一個細胞陣列



產生細胞陣列可利用指定陳述法,或用cell指令先安排細胞陣列之大小,然後再將資料移轉至各細胞存放。指定陳述法係將資料直接指定給特定細胞元素,一次一個。 MATLAB 將自動建立對應之細胞陣列。

將資料指定給某特定細胞有兩種方式:即細胞索引(cell indexing)與內容索引(content indexing)。細胞索引是先在等號左邊以下標說明細胞之位置,並且用一般標準的括符()括起來,如A(2,3)是;而等號之右邊則必須將所要設定的內容利用大括等{}括起來,表示該位置所設定的內容。例如為產生一個2 x 2的細胞陣列A:

>>A(1,1) = {[3 4 5; 1 5 7; 3 5 4]};
>>A(1,2) = {'馮丁樹'};
>>A(2,1) = {2:3:11};
>>A(2,2) = {-5+4j};
>>A
A =
[3x3 double] '馮丁樹'
[1x4 double] [-5.0000+ 4.0000i]

由這裡可以看出,產生一個細胞陣列可以如一般的矩陣指定方式一樣,將內容逐一填入。由於細胞陣列A之內容格式可能不一,故任何資料均可存放。此外,細胞中若為矩陣,則僅以矩陣之大小及數值內容表示。要注意的是等號右邊需用大括符,不能使用矩陣的中括符。而{}表示為空集合,是一種合法的型式。若要知道其實際內容,可以使用celldisp指令:

>>celldisp(A)
A{1,1} =
3 4 5
1 5 7
3 5 4

A{2,1} =
2 5 8 11
A{1,2} =馮丁樹
A{2,2} = -5.0000 + 4.0000i

此時所有內容均可以列出。若要以圖示表示此細胞陣列之相關位置,亦可使用cellplot指令,這是一個繪圖式指令。它會將各細胞之元素用圖表示出,但矩陣之內容仍不顯示。

>>cellplot(A)




內容索引則是在等號左邊以大括符定其位置,在等號右邊則置入實際之內容。其方式如下:

>>A{1,1} = [3 4 5; 1 5 7; 3 5 4];
A{1,2} = '馮丁樹';
A{2,1} = [2:3:11];
A{2,2} = -5+4j;
>>A
A =
[3x3 double] '馮丁樹'
[1x4 double] [-5.0000+ 4.0000i]

結果相同。換言之,在設定時,大括符僅能使用一次,不是在等號之左邊,就是在右邊。若所設定的資料位址超過原來的範圍,該細胞陣列會自動擴張,以配合需求。例如前述之A細胞陣列大小為2 x 2,若另外加一項內容,則會變為 3x3:

>>A(3,3)={'美麗的夏天'}
A =
[3x3 double] '馮丁樹' []
[1x4 double] [-5.0000+ 4.0000i] []
[] [] '美麗的夏天'

細胞陣列的語法中,利用大括符形成細胞內容之結構,與中括符作為一般矩陣的結構功能一樣;唯一相異之處是大括符允許其內又有大括符的巢狀結構。因此,細胞陣列中亦可包括另一個細胞陣列。在大括符中,係利用逗點或空格表示行向斷點,用分號作為細胞間之列向斷點。例如:

>>C = {[5 1 2], [9 7]; [1 7 3 4], [8 6 3]}
C =
[1x3 double] [1x2 double]
[1x4 double] [1x3 double]

利用cell指令則可預先設定一個細胞陣列之空集合,以事先保留此細胞陣列之架構。其中cell(n)可產生一個n x n之空細胞陣列:

>>D=cell(3, 4)
D =
[] [] [] []
[] [] [] []
[] [] [] []

若採用多項輸入時,可以直接利用陣列的內容陳述,元素中若為矩陣則必須以矩陣之方式表示;若為文字,則必須加上撇號,其前後則用大括符括起來:

>>D={'天龍八部' '金庸' 1988 [200 220 450]
'紅樓夢' '曹雪芹' 1890 [150 180 160]
'西遊記' '吳承恩' 1880 [120 140 100]};

>>D
D =
'天龍八部' '金庸' [1988] [1x3 double]
'紅樓夢' '曹雪芹' [1890] [1x3 double]
'西遊記' '吳承恩' [1880] [1x3 double]

若要單項輸入或修改,則在等號左邊採用一般括符之索引;右邊則採用大括符置放每個細胞之資料:

>>D(3,3)={1885}
D =
'天龍八部' '金庸' [1988] [1x3 double]
'紅樓夢' '曹雪芹' [1890] [1x3 double]
'西遊記' '吳承恩' [1885] [1x3 double]

實際上其指定方式亦可依細胞陣列元素之總排序之行列順序(亦即依第一行、第二行、…等之順序)。就上述第(3,3)之位置而言,其總排序順位應為9,故:

>>D(9)={1880}
D =
'天龍八部' '金庸' [1988] [1x3 double]
'紅樓夢' '曹雪芹' [1890] [1x3 double]
'西遊記' '吳承恩' [1880] [1x3 double]

又可將西遊記的年代更正為1880了。若要更改某細胞元素內之矩陣中之特定元素,則必須先使用大括符{}先指定細胞之位址,再用標準括符()定矩陣之位置。例如要更改D(3,4)中的第二項140為240時:

>>D{3,4}(1,2)=240
D{3,4}
ans =
120 240 100

顯示該項已經改正。

呼叫的程序亦相同,例如位置D(2,2):

>>D(2,2)
ans =
'曹雪芹'

>>D(1:3)
ans =
'天龍八部' '紅樓夢' '西遊記'

>>D(7:9)
ans =
[1988] [1890] [1880]

注意使用標準括號僅取其位址,而非其內容,要求第9項之內容必須使用大括符:

>>D{9}
ans = 1880

使用大括符可以將其內容組成新的矩陣,以供後來之應用:

>>P=[D{10};D{11};D{12}]
P =
200 220 450
150 180 160
120 240 100

利用細胞索引可以取出次陣列,設定給另外變數:

>>B=D(:,1:3)
B =
'天龍八部' '金庸' [1988]
'紅樓夢' '曹雪芹' [1890]
'西遊記' '吳承恩' [1880]

若僅為單細胞元素之內容,則使用一般括號或大括號索引均可,但若要取其中之之一項或部份時,則必須先使用大括符定位,然後再用一般型式索取矩陣中之個別內容。

若要取得細胞陣列中某特定項之內容,並置放對應之新變數時,可使用deal指令,但兩邊之項目必須相對應:

>>[book, author]=deal(D(:,1),D(:,2))
book =
'天龍八部'
'紅樓夢'
'西遊記'

author =
'金庸'
'曹雪芹'
'吳承恩'


deal這個指令等於進行配對的意思,將等號右邊之對應的資料配給左邊的變數,有如等號,只是其功能在於可同時處理許多變數或一系列變數,故特別適用於細胞陣列或下節要談到的結構資料。若右邊的資料僅有一項,則會以同樣的資料配給左邊的所有參數。這些參數可能是獨立的參數,也可能是另一陣列。例如:


>> sys = {rand(3) ones(3,1) eye(3) zeros(3,1)};
>> sys

sys =

[3x3 double] [3x1 double] [3x3 double] [3x1 double]

>> [a,b,c,d] = deal(sys{:});
>> a

a =

0.9501 0.4860 0.4565
0.2311 0.8913 0.0185
0.6068 0.7621 0.8214

>> b

b =

1
1
1

>> c

c =

1 0 0
0 1 0
0 0 1

>> d

d =

0
0
0



在細胞陣列中,若要刪除某一細胞,只要送給該細胞一個空集合即可,即D(…)=[]。這種細胞陣列亦可利用reshape函數進行大小重排,但整排後之矩陣大小需為原矩陣可除盡之數。

實際上細胞陣列有什麼好處?基於其容許不同大小及型式的內容之特性,有些應用依需要而定。它可以取代一些以逗點分開的變數名單,其中包括函數輸出入名單、作業過程之展示及陣列之建構等。由於MATLAB處理細胞陣列有如分離的變數。如:

>>C={[1 2 1], [1 0 1]}
C =
[1x3 double] [1x3 double]

可以使用細胞陣列中之元素作為其他函數之輸入參數,例如將兩多項式之係數相乘:

>>d=conv(C{1},C{2})
d =
1 2 2 2 1

5.4 結構陣列

5.4 Structural Array


MATLAB之第五版以上也具有結構型之資料型式,可以利用欄位名稱輸入資料庫或直接呼叫,對於以後寫資料庫之擷取型式具有相當大的方便性。而結構陣列之存在,亦可不拘資料之屬性,可為一搬文字敘述或矩陣之型式。這種結構與一般資料庫之分類法則相同,應用時相當方便。


patient
.name ‘張三豐’
.billing 3,838.00
.test [12 01 70 180
3 20 65 150
5 30 79 150]



結構陣列之型式以一個病歷作說明:在一個資料庫下有欄位名稱,若以patient為資料庫之總稱時,其下有name、billing及test等三個欄位,第一個欄位為病人之姓名,第二欄為帳單,第三欄為其檢驗之記錄;後者因每一次有相同的實驗項目,故以矩陣表示。

要建立這些資料並不困難,可在指令窗下輸入:


>>patient.name='張三豐';
>>patient.billing=3838.00;
>>patient.test=[12 01 70 180;3 20 65 150;5 30 79 150];


輸入之後,只要打入patient,即有這筆記錄:


>>patient
patient =
name: '張三豐'
billing: 3838
test: [3x4 double]


最後一項是一矩陣,它僅顯示該矩陣之大小,未顯示其內容,若要知道其內容,可以下指令:


>>patient.test
ans =
12 1 70 180
3 20 65 150
5 30 79 150


由此可見,要知道其下之欄位內容,可以在名稱與欄位名稱間加一點即可,如:


>>patient.name
ans =張三豐
>>patient.billing
ans = 3838


在其他運算式中,若要參照某一欄位內容之值時,亦可採用同樣的語法取得該值。

上面之結構仍然僅能指向某一位病人而已,若有其他病人也需同樣的結構如何表示呢?此時必須在病人之參數加上指標,以表示不同的病人利用同樣的結構建檔。如下建立第二筆記錄:


>>patient(2).name='李四文';
>>patient(2).billing=2424.00;
>>patient(2).test=[11 05 75 200;4 20 80 160;6 30 55 140];


由於patient之內容不只一個,此時若要知道此一patient的資料如何,仍然可在指令窗下指令:


>>patient
patient =
1x2 struct array with fields:
name
billing
test


但其結果出人意外的與前面不同,主要是我們已將patient變成一個1x2 的結構性矩陣,而其下之欄位名稱則會加以標明,不再標出其內容,不過若下如下之指令,則會有所要的資料:


>>patient(2)
ans =
name: '李四文'
billing: 2424
test: [3x4 double]

>>patient(1)
ans =
name: '張三豐'
billing: 3838
test: [3x4 double]

至少目前應有張三豐及李四文兩筆病人的記錄。其餘可以採用類似方法增加筆數。基本上每一新記錄均需維持相同的欄位名稱及數量,才能組成一個結構。

雖然這些資料均以結構型式存於檔案中,但如前面所述,其單獨之變數仍可以取出隨意運用,如此可以增加變數應用之彈性,例如繪出張三豐之實驗之資料:

>>bar(patient(1).test)

張三豐的病歷資料



故其內在之資料要如何處理,其情況與一般的手法相同,只是目前仍然不知道張三豐的病歷到底表示的是什麼而已。

若要瞭解這個檔案之欄位名稱,除直接打入該檔的名稱外,亦可用fieldnames這一個指令:


>>fieldnames(patient)
ans =
'name'
'billing'
'test'


其中之某欄位亦可利用rmfield這個指令去除掉:


>>patient2=rmfield(patient,'test')
patient2 =
1x2 struct array with fields:
name
billing


因此,基本結構patient2與patient相同,但少了test這個欄位。一次若要除去好幾個欄位時,則必須使用文字矩陣來處理。以往討論到矩陣時,均以數值為主,故矩陣內的項目均為數值格式,且用中括號[]作為範圍;碰到項目為文字型式時,其處理方式需以單一之文字為對象。若要採用長度不一的文字矩陣時,則各項之文字必須用撇號 ’ 號括起來,其矩陣亦必須以大括號{}作為範圍,如{'張三豐','李四文','陳五分'}是。故若要除去許多欄位時,可採用下面的方法:


>>patient2=rmfield(patient,{'test' 'billing'})
patient2 =
1x2 struct array with fields:
name


因此,patient2僅剩下name一欄了。當然欄位亦可無中生有,只要將要直接添加新欄位並設定一個值給它即可,如:


>>patient2(1).billing=20
patient2 =
1x2 struct array with fields:
name
billing


換言之,整個檔案均會增加一個欄位billing,而第二筆記錄欄位已存在,但仍為空值,因為尚有未該筆記錄之billing的內容:


>>patient2(2).billing
ans =
[]

>>patient2(1).billing
ans =
20


上述方法亦可用另一指令setfield來逹成:


>>setfield(patient2(2),'billing', 30)
ans =
name: '李四文'
billing: 30


事實上,結構化檔案亦可以利用struct指令來達成:


>>patient3=struct('name',{'陳五分','趙六錢'},'billing',{5656.00, 7878.00},'test',{[1 2],[3 4]})
patient3 =
1x2 struct array with fields:
name
billing
test

>>patient3(1)
ans =
name: '陳五分'
billing: 5656
test: [1 2]

>>patient3(2)
ans =
name: '趙六錢'
billing: 7878
test: [3 4]


上述指令之結構可能需要補充說明一下,以便發揮一指神功的效果。struct後面括號內包含的項目為各欄位及對應之值。當然這些值可能為文字,也可能為數值。基本上文字均需使用撇號 ’ 括起來。對應欄位後面應使用中括號[]括起來,裡面的項目內容若有好幾筆資料,則需用逗號分開。各筆若為矩陣,則需依矩陣的表示法輸入,並用中括號[]括起來。下面兩個例子可能給大家一些概念:


>>s = struct('greetings',{{'Hello, there','新年好,萬事如意'}},'lengths',[5 3])
s =
greetings: {'Hello, there' '新年好,萬事如意'}
lengths: [5 3]

>>s(1)
ans =
greetings: {'Hello, there' '新年好,萬事如意'}
lengths: [5 3]


顯然s內容僅有一項,但greetings內應有兩項:


>>s.greetings(1)
ans =
'Hello, there'

>>s.greetings(2)
ans =
'新年好,萬事如意'


再看看同一個例子,但稍做修改:


>>s = struct('greetings',{'Hello, there','新年好,萬事如意 '},'lengths',[5 3])
s =
1x2 struct array with fields:
greetings
lengths

>>s(1)
ans =
greetings: 'Hello, there'
lengths: [5 3]

>>s(2)
ans =
greetings: '新年好,萬事如意 '
lengths: [5 3]


此時s有兩筆記錄,且其lengths:均為 [5 3],足見不足的項目它會用同一個值補上去。
再看看它在說明中提供的例子:


>>ss = struct('type',{'big','little'},'color','red','x',{3 4})
ss =
1x2 struct array with fields:
type
color
x

>>ss(1)
ans =
type: 'big'
color: 'red'
x: 3

>>ss(2)
ans =
type: 'little'
color: 'red'
x: 4


這樣的結果會更清楚了,讀者應可從中體會其功能。下面是一個比較複雜的例子:


>>A=struct('data',{[3 4 7;8 0 1], [9 3 2; 7 6 5]}, ...
'personal',{struct('testyear',1999,'height',[163 178 172],...
'weight',[80 86 60]),struct('testyear',2000,'height',[157 168 159],...
'weight',[48 50 48])})
A =
1x2 struct array with fields:
data
personal

>>A(1)
ans =
data: [2x3 double]
personal: [1x1 struct]

>>A(1).data
ans =
3 4 7
8 0 1

>>A(1).personal
ans =
testyear: 1999
height: [163 178 172]
weight: [80 86 60]

>>A(2).personal
ans =
testyear: 2000
height: [157 168 159]
weight: [48 50 48]


其他相關指令如disp也可以看其結構:



>>disp(A)
1x2 struct array with fields:
data
personal


而函數isfield則可檢驗檔案結構內是否有某一特定欄位,如:


>>f=isfield(A,'data')
f =
1

下表為有關處理matlab結構之各函數:







dealDeal inputs to outputs.
fieldnamesGet structure field names.
isfieldReturn true if the field is in a structure array.
isstructReturn true for structures.
rmfieldRemove a structure field.
structCreate or convert to a structure array.
struct2cell Convert a structure array into a cell array.

10/20/2006

第四章 內建函數

函數是針對某特定功能所撰寫的程式,利用其固定的名稱執行所需之任務。函數通常可分為內建函數與自訂函數兩大類。內建函數則又視其公用程度而分一般函數與特定函數。目前我們討論者,以一般函數為原則。茲分別作說明:

4.1數學基本函數

一般函數有三角函數、矩陣函數及統計函數。有些由於解說之需要,在前文中已有述及。但基本上,其歸類有些也甚為模糊的,只是這種分類也沒多大意義。故讀者參考時,可以前後作對照。三角函數是工程上最為常見的函數,下表將其列入。

4.1.1三角函數


由於三角函數牽涉到角度的單位,以往其輸入參數均以弧度為單位,故應用時須先轉換,令pi等於180度,依此比例計算。新版的應用中,則增加以度數為單位輸出與輸入的型式,因此不必再經過轉換的程序。在函數的名稱上,若輸入為度數,則可使用字尾加"d"的型式,如xxxd是,詳情見表4.1。當然是否必須使用不同的度數為指令,則見仁見智。若考慮新舊版的共通性,仍仍以採用轉換之型式為宜。

表4.1三角函數指令表













指令型式說明
sin / sind正弦函數,輸入參數單位為弧度/度數
cos /cosd 餘弦函數,輸入參數單位為弧度/度數
tan /tand 正切函數,輸入參數單位為弧度/度數
cot /cotd 餘切函數,輸入參數單位為弧度/度數
asin / asind 反正弦函數,輸出參數單位為弧度/度數
acos /acosd 反餘弦函數,輸出參數單位為弧度/度數
atan /atand 反正切函數,輸出參數單位為弧度/度數
acot /acotd 反餘切函數,輸出參數單位為弧度/度數
sinh / cosh 超正/餘弦函數
asinh /acosh 反超正/餘弦函數,輸出參數單位為弧度
tanh /coth 超正/餘切函數
atanh /acoth 反超正/餘切函數,輸出參數單位為弧度
atan2 四象限反正切函數,輸出參數單位為弧度


例1. 計算sin(60)之值:

>>sin(60*pi/180)

ans = 0.8660

>>sind(60)

ans = 0.8660



其結果相同。

例2. 計算sin²(60度)+cos²(60度)之值

>>theta=60*pi/180;
>>sin(theta)^2;+cos(theta)^2;
ans = 1

>>sind(60)^2;+cosd(60)^2;
ans = 1


兩者之結果亦相同。

例3. 計算arctan(1)+arccot(1)之值

>>atan(1)+acot(1)
ans = 1.5708

>>atand(1)+acotd(1)
ans = 90


前者之單位為弧度,後者為度數,其實兩者所指的角度是相同的。

表4.1中亦列有超函數之計算,這是在許多工程數學常見到的解,其型式近乎雙曲線函數。其定義如下:

cosh x ={exp(x)+exp(-x)}/2

sinh x ={exp(x)-exp(-x)}/2


若依其定義進行計算,則可利用後節之對數函數配合解決,若直接用超函數表示,則可簡單計算如下:

>>cosh(1)
ans =
1.5431


4.1.2指數與對數函數



有關指數與對數函數之指令則如表4.2所示。這些指令與極座標有關,也與超函數之運算有關。

表4.2 指數與對數與複數函數指令表











指令型式說明
sqrt(x)開方根( )
exp(x)指數函數(ex)
log(x)自然對數( )
log10(x)以10為基底之對數(log10x)
abs(x)(複數之)絕對值
angle(x)複數或極座標之夾角,弧度
real(x)複數之實數部份
imag(x)複數之虛數部份
conj(x)共軛複數
i,j虛數單位,如3+4j,或3+4i


一般數之次方常用"^"表示,如A^5代表A5。開平方根則是常用的指令,故有特別的操作函數sqrt,例如3的開平方:

>>sqrt(3)
ans = 1.7321


同理負數的開方也可以得到結果,只是會產生複數,以-5之開方為例:

>>sqrt(-5)
ans = 0 + 2.2361i


結果進入虛數的世界,因為得到虛根。所以MATLAB在處理數值時,開方根並不限於正值。
指數方面,因為基底之不同,數學上分為以e 與以10為底;對應於指令方面,分別為log(x)與log10(x)之型式。與對數相應,以e 為底者稱為自然對數,用exp(x)表示。例如e¹其值為:


>>exp(1)
ans =
2.7183


顯然e值就是2.7183。由於定義上,y=e^x 兩邊分取對數,則log(y)=x。故上述結果代入log(y)應可得到x,如:

>>log(2.7183)
ans =
1.0000


所以兩者是相互為用的。在以10為基底的情況下,若要更換為自然對數,則其對數值應乘以log(10) 。以x=10^5為例,兩邊取自然對數後,成為log(x)=5*log(10),則:


>>logx=5*log(10)
logx = 11.5129

>>x=exp(logx)
x =
1.0000e+005


其值與原數相同。注意:結果為以10為底之指數型式,由e後面為其次方值。
指數的型式,其輸入項亦可使用複數,以配合極座標之應用。通常將複數中之實數部份作為橫軸,虛數之值則視為縱軸,如此可以構成二度空間。假設某向量R=3+4j,則:


>>R=3+4j
R = 3.0000 + 4.0000i


由此可以利用計算其絕對值abs(R)、夾角angle(R):

>>abs(R)
ans = 5

>>angle(R)*180/pi
ans = 53.1301


後者之角度已轉換為度數。利用同一複數量R亦可求其實數、虛數及共軛複數:

>>real(R) %實數部份 
ans = 3

>>imag(R) %虛數部份
ans = 4

>>conj(R) %共軛複數
ans = 3.0000 - 4.0000i


根據尤拉公式,虛數之自然對數可以與正、餘弦作轉換,其基本型式為:

eix=cos(x)+ i sin(x)

此時,x之值為弧度值。因此利用這種轉換亦可將自然對數之型式轉換為三函數之關係:例如,設theta=[ 20 60 80]度,則

>>theta=[20 60 80]/180*pi %先轉換為弧度
theta = 0.3491 1.0472 1.3963

>>eix=exp(i*theta)
eix = 0.9397 + 0.3420i 0.5000 + 0.8660i 0.1736 + 0.9848i

>>cos(theta)+i*sin(theta)
ans = 0.9397 + 0.3420i 0.5000 + 0.8660i 0.1736 + 0.9848i


由此可以證明三個不同角度下,其演算結果相同。

就超函數之定義而言,其值與指數間具確定的關係,依照前述之cosh(x)與sinh(x)之定義,兩者相加可以得到如下之如果:


ex=cosh(x)+sinh(x)

>>cosh(1:5)+sinh(1:5)
ans = 2.7183 7.3891 20.0855 54.5982 148.4132

>>exp(1:5)
ans = 2.7183 7.3891 20.0855 54.5982 148.4132



兩者對應之結果一致。

4.1.3元素操作函數


有些數值必須在計算中進行整理,使其符合需求。這些函數包括四捨五入(round)、頂部整數(ceil)、底部整數(fix)、近零整數(floor)等。這些功能如下表:

表4.3元素操作函數






指令型式說明
ceil(x)大於x之最近整數
fix(x)位於接近於零之整數
floor(x)小於x之最近整數
round(x)四捨五入之整數
sign(x)符號值,若x>0為+1;x=0為0;x<0為-1>


假設有一向量A為:


>>A=rand(1,10)*10-5 %利用隨機函數產生-5至+5間之亂數
A =
1.1543 2.9194 4.2181 2.3821 -3.2373 -0.9429
   4.3547 4.1690 -0.8973 3.9365

>>round(A) %進行四捨五入。
ans =
1 3 4 2 -3 -1 4 4 -1 4

>>fix(A) %最接近零之整數
ans =
1 2 4 2 -3 0 4 4 0 3

>>floor(A) %得到最接近A值但低於A值之整數
ans =
1 2 4 2 -4 -1 4 4 -1 3

>>ceil(A) %得到大於A且最接近於A之整數
ans =
2 3 5 3 -3 0 5 5 0 4

>>sign(A) %得到該值之符號值,若該值為零,則為零(本例沒有)
ans =
1 1 1 1 -1 -1 1 1 -1 1

4.2矩陣之統計

對於一個矩陣元素值之統計,前面也有部份討論到,諸如sum、mean、max、min等等,都是就元素值之特性進行操作者。其中有些以向量矩陣為主(包括行向量與列向量),但亦可應用於一般矩陣(例如m x n矩陣),將每行之元素進行處理(例如累加),然後將結果置於一個列矩陣中。這些統計型函數如表4.5所示:

表4.5矩陣簡單之統計指令








指令型式說明
median(A)求陣列A之中間值
cumsum(A)求陣列A之累加值
prod(A)求陣列A之乘積值
cumprod(A)求陣列A之累積值
std(A)求陣列A之標準差




4.2.1 MEDIAN中間函數


求矩陣中之中間值數,其與mean函數不同,後者求其平均值,但應用上大致相同。


>>A=magic(5)
A =
17 24 1 8 15
23 5 7 14 16
4 6 13 20 22
10 12 19 21 3
11 18 25 2 9

>>median(A)
ans =
11 12 13 14 15


此為各行之中間值,若要求各列之中間值則可用median(A,2)



>>median(A,2)
ans =
15
14
13
12
11



4.2.2 CUMSUM累加函數


這是矩陣之累加函數,可以就行向各元素逐一累加。若仍以上式A魔術矩陣為例:


>>cumsum(A)
ans =
17 24 1 8 15
40 29 8 22 31
44 35 21 42 53
54 47 40 63 56
65 65 65 65 65


其所得結果為由上往下累加,若要列向累加,則可以在其後加一個dim=2之參數:



>>cumsum(A, 2)
ans =
17 41 42 50 65
23 28 35 49 65
4 10 23 43 65
10 22 41 62 65
11 29 54 56 65


而魔術方陣其各方向累加起來之值均相同,也可由此得證。利用cumsum函數也可以立即產生一個等差級數之陣列,例如等差值5之陣列如下:



>>M=cumsum(0:5:50)
M =
0 5 15 30 50 75 105 140 180 225 275




4.2.3 PROD乘積函數


prod函數與sum之功能大體相同,只是prod是將各元素沿特定陣列方向相乘,得其乘積為答案。例如:



>>prod(1:5)
ans = 120


將上式之魔術方陣A進行行向運算:



>>prod(A)
ans =
172040 155520 43225 94080 142560


顯然,其各行結果不像加法,會有相同的結果。若要列向相乘,也可以設dim=2:



>>prod(A,2)
ans =
48960
180320
137280
143640
89100



4.2.4 CUMPROD累積函數


累積函數cumprod與累加 cumsum亦屬相同的功能,只是以乘法表示而已。其可能應用的則為等比級數之產生。例如:



>>cumprod(1:2:10)
ans =
1 3 15 105 945



4.2.5 STD標準差


一般標準差之計算有兩種方式,其一為分母除以n-1,即不包括平均值之項目,一為分母除以n,
計算上式標準差之函數型式如下:







s=std(X, flag, dim)


其中dim之定義與前面相同外,flag之旗標為0則採用第一式計算標準差,換言之,其分母為n-1;若旗標為1,即採用第二式計算。一般採用的方式以旗標為0者為多,其效果與不標旗標的情況相同。相關例子如下:設M為隨機函數:



>>M=round(rand(4,5)*100)-50
M =
34 33 -20 -20 -12
-48 0 -31 4 36
18 21 -31 -35 35
-12 -7 18 20 9

>>std(M)
ans =
36.1109 18.5000 23.2522 24.5000 23.0217

>>std(M,1)
ans =
31.2730 16.0215 20.1370 21.2176 19.9374


std(M)與std(M,1)之差異在於前者除以n-1,後者除以n。在項目較少的場合,後者之變異數會較小,在統計學上會有顯著的差異。通常標準差常與平均值相配合,其平均值可以對應求得如下:



>>mean(M)
ans =
-2.0000 11.7500 -16.0000 -7.7500 17.0000


若要求得列向的標準差,則可令第三個參數dim=2,得到的是行向量標準差。



>>std(M,0,2)
ans =
28.0357
32.6833
32.2614
14.5017


若需求得整個矩陣之平均值與標準差,則可改用下面方式計算:



>>[mean(M(:)) std(M(:))]
ans =
0.6000 26.0897


由於M矩陣為隨機函數產生,故其平均值應接近於零。

4.3矩陣之應用

矩陣之相乘雖然屬於元素與元素間之相乘,但其結果可能因解釋或其實際之意義而不同。前文所述之陣列相乘即是一種例子,但其結果是相對應位置的元素相乘。這種逐元乘法在Matlab中相當普遍,而通常矩陣的加法、減法是與逐元之加法、減法是一致的,其應用比較簡單。至於矩陣的乘法(或未來的除法)其代表的意義自有不同。

範例一:向量之點積(Dot product)


以向量為例,通常一個向量可以使用三個垂向座標來表示,而各座標之單位向量(其量為一,但方向維持座標方向)若以i、j、k表示,則任一向量V、U可以表示如下:

V=v1 i + v2 j + v3 k;U=u1 i + u2 j + u3 k


今若要求這兩向量間之內積(dot product),則V。U之結果應為各向量同方向之分量之積之和,亦即:


V。U=v1.u1+v2.u2+v3.u3

這個結果若能以矩陣表示,設V=[ v1 v2 v3];U=[ u1 u2 u3] ,則

V。U=v1.u1+v2.u2+v3.u3= [ v1 v2 v3][ u1 u2 u3]'

在MATLAB有一個點積指令稱為dot(),其指令型式如下:

dot(U,V)

假設:V=[ 25 4 50]; U=[15 20 5],結果:


>>V.U=[25 4 50]*[15 20 5]'
ans =
705
>> [15 20 5]*[25 4 50]' %即使順序掉換,其結果亦相同。
ans =
705

>> [25 4 50]'*[15 20 5] %但是這樣掉換,其結果完全不同。
ans =
375 500 125
60 80 20
750 1000 250


上述的例子中,使用點積指令亦可得到同樣的結果:


>> dot([25 4 50],[15 20 5])

ans =

705

>> dot([25 4 50],[15 20 5]')

ans =

705



因此,若利用點積指令,也可以求得某特定向量之絕對值如下:


>> r=sqrt(dot([25 4 50],[25 4 50]))

r =

56.0446

>> u=[25 4 50]/r

u =

0.4461 0.0714 0.8921


上述u則為該向量之單位向量了。

範例二:叉積指令(Cross)


兩向量間之叉積在於求得兩向量所構成之垂向面積,而其方向則是垂直於該面積之正方向。例如扭矩等於力與力臂之叉積,而其方向則在於共同垂直於此二向量之方向。叉積之指令型式為:

M=cross(R,F)

設有A向量沿x軸,其向量值為A=1i,以向量陣列表示為A=[1 0 0];B為沿y方向,其向量值B=2j,或B=[0 2 0];則結果向量D應為:

D=AxB (注意此處x號為叉積或cross的意思,不要與矩陣之*混淆),則


>> A=[1 0 0];B=[0 2 0];
>> D=cross(A,B)

D =

0 0 2


其結果在z方向,其值為2。
若兩向量均不在特別之座標上,則結果也可能在空間上,例如:


>> m=cross([1 2 3],[5 4 3])

ans =

-6 12 -6


至於要求得m之絕對值,則有下列方法:


>> r2=sqrt(m*m')

r2 =

14.6969

>> r2=sqrt(dot(m,m))

r2 =

14.6969


兩者均可得到相同的答案,則如何求其單位向量,相信應更容易了。

範例三:九九乘法表


例如:製造一個九九乘法表,其方法可以利用上述的用法迅速得到答案:

>> A=[1:9];Table=A'*A
Table =
1 2 3 4 5 6 7 8 9
2 4 6 8 10 12 14 16 18
3 6 9 12 15 18 21 24 27
4 8 12 16 20 24 28 32 36
5 10 15 20 25 30 35 40 45
6 12 18 24 30 36 42 48 54
7 14 21 28 35 42 49 56 63
8 16 24 32 40 48 56 64 72
9 18 27 36 45 54 63 72 81



同理,試算出一長L=[2 4 6 8 10]m、寬W=[1 2 3 4 5]m所構成之組合及其對應之重量。設密度為10kg/cm^3 。


>> Weight=[2:2:8]'*[1:5]*10

Weight =
20 40 60 80 100
40 80 120 160 200
60 120 180 240 300
80 160 240 320 400


範例四:生產成本分析


某工廠件生產三件主力產品A、B、C三種。生產過程中,其材料、勞力、設備及運輸等之單位產品成本亦自不同。設每樣產品一年分四季生產,其售價亦四季不同,求其每季之生產成本、單項產品之生產成本、單項產品之銷售及利潤等。其相關資料如下:

單位生產成本,P(NT$)




產品名材料 設備人力運輸
100226512
5015805
200305020


每季生產量,Q單位




產品名第一季第二季第三季第四季
10020025090
250200150100
80100120150


每季售價格,W(NT$)




產品名第一季第二季第三季第四季
300250250200
400350250150
500600650700


先將三種資料輸入P、Q、W等三個矩陣,其大小均為(3x4):但Q之四項與P之四項並不屬同性質,故即使矩陣之大小配合,也須注意不能隨意相乘。


>>P=[100 22 65 12;50 15 80 5;200 30 50 20];
>>Q=[100 200 250 90;250 200 150 100;80 100 120 150];
>>W=[300 250 250 200; 400 350 250 150;500 600 650 700];


輸入後之值如下:


>>P %單位成本
P =
100 22 65 12
50 15 80 5
200 30 50 20 Q %每季生產量

>>Q =
100 200 250 90
250 200 150 100
80 100 120 150 W %每季之產品價格

>>W =
300 250 250 200
400 350 250 150
500 600 650 700



為求得每季之生產成本,可用C=P'*Q 計算,其大小之型式為(4x4)=(4x3)(3x4)。結果為4x4,其列為材料、設備、人力及運輸之成本,行向為四季。注意前者P用移置功能,主要將其大小轉為符合矩陣之乘法:


>>C=P'*Q
C =
38500 50000 56500 44000
8350 10400 11350 7980
30500 34000 34250 21350
4050 5400 6150 4580



上式C矩陣為四季之成本,若將材料、設備、人力及運輸成本合計,則各季之成本小計可用sum函數計算:


>>Qcost=sum(C)
Qcost =
81400 99800 108250 77910


Qcost為四季之成本,若需要總成本,則需再用一次sum函數:


>>Total_cost=sum(Qcost)
Total_cost =
367360


其次要計算單項產品之成本時,可利用原來之P矩陣。因為其列為產品名稱,行為分項成于,故若列向相加,應可得各項產品之單位成本,其法可利用sum函數,但其第二參數應輸入2,表示以列向(橫向)相加:


>>Prod_unitcost=sum(P,2)
Prod_unitcost =
199
150
300


同樣,對生產量Q而言,由於其行向屬不同季節,若將其合併,應為一年之總生產量,故仍然應用sum函數,但記得其第二參數為2:


>>Prod_Q=sum(Q,2)
>>Prod_Q =
640
700
450


有了單位總生產成本Prod_unitcost及年生產量Prod_Q,即可直接用元素對應元素乘法,求得其總成本。注意此時之操作元乘號前須加點,即.* 。


>>Prod_cost=Prod_unitcost.*Prod_Q
Prod_cost =
127360
105000
135000


Prod_cost 即為年生產總成本。其次要求得年銷售金額,此時因為Q與W矩陣之行向同為季別,故只要使用元素對元素之乘法即可。其四季之銷售量D=Q.*W:


>>D=Q.*W
D =
30000 50000 62500 18000
100000 70000 37500 15000
40000 60000 78000 105000


目前為求得各項產品之銷售金額,故只要每項產品之四季金額相加即可,此時仍使用sum函數,其第二輸入參數仍設定為2:


>>Prod_sale=sum(D,2)
Prod_sale =
160500
222500
283000


Prod_sale為全年三種產品之銷售金額。扣除上項得到之成本Prod_cost,即為利潤Profit:


>>Profit=Prod_sale-Prod_cost
Profit =
33140
117500
148000


綜合判斷,C產品之獲利較高;B產品次之。讀者可就四季所得之利潤進行分析,其結果應更具參考性。

10/09/2006

第三章 程式結構與設計

前面介紹指令及函數檔的過程中,已陸續介紹一些程式指令之用法及功能。基本上程式是一些指令之集合,執行後可以獲得所希望的結果。複雜的程式則是採用主程式、副程式或函數檔間相互呼叫的方式,由各檔案專司解決特定的功能,如此形成一個程式之大架構。這些程式執行過程中,必須配合執行之中間值,進行比較、判斷,最後得到適當之結果。有些檔案可能需經過多次呼叫,才能達到程式預設之目的;有些則需判斷結果之適當性,依據最佳的決策選擇執行的路徑。

在程式結構中,一般除基本指令外,尚須配合流程之控制指令,使程式能配合特定條件運作。控制指令以邏輯值進行比較,其操作元包括關係操作元及邏輯操作元等,可以應用於變數及矩陣。在程式之運算過程中,為執行某一特定任務均有其基本之演算法。其中必須藉助控制指令以改變指令執行之順序,或稱為控制結構。這些結構必須藉其特定之關係比較,就比較結果決定程式進行的順序。這些運算法大體上可分為循序運算、條件運算及迴圈運算等三種方式。本節將針對這三種特定型式之指令作詳細的介紹。

3.1關係操作元

兩陣列間之相互比較是常有的事,若能配合程式之關係運算則可節省撰寫時間。一般關係操作元主要在兩個變數間作邏輯性比較,,然後選擇一個最佳的結果;在MATLAB之應用上實也脫離不了這些原則。只是在陣列之比較上,必須時時注意其大小之匹配。下面為一般操作元之定義:


表3.1 關係操作元










操作單元符號說明
<小於
<=小於或等於
>大於
>=大於或等於
==等於
~=不等於



上述關係操作元當中,比較特別需要注意的是,相等之關係操作元必須使用兩個等於,以便與設定變數值之等於有所區別。

在同樣大小的矩陣中,其比較可依元素與元素相互比較,將結果置於另一矩陣內。以下面兩矩陣A、B為例,兩者分別由魔方陣及亂數函數產生,大小均為4x4:

>>A=magic(4)
A =
16 2 3 13
5 11 10 8
9 7 6 12
4 14 15 1

>>B=round(rand(4)*16)
B =
7 14 13 13
15 8 0 8
7 3 11 11
7 11 6 7

>>C=A>=B
C =
1 0 0 1
0 1 1 1
1 1 0 1
0 1 1 0

新的矩陣C為兩個矩陣比較後之結果,其元素由邏輯值0或1組成,大於等於為1否則為0。

在矩陣之操作中,條件型態也可以加入,以歸納一些符合條件的元素。以上述A矩陣為例,若在A矩陣中,將所有A>=B之元素均設為100,則可利用上述C矩陣之邏輯值作為條件如下:

>>A(C)=100
A =
100 2 3 100
5 100 100 100
100 100 6 100
4 100 100 1

A矩陣中只要相對應於C矩陣中之值為真,或合於A>=B條件之元素因而均被設為100。這樣的結果僅需一個指令即可達成。實際上要達成上項目的並不一定要經過C矩陣,直接採用下列方式亦可:

>>A(A>=B)=100
>>A =
100 2 3 100
5 100 100 100
100 100 6 100
4 100 100 1

其結果與上述方法相同。這種方法並不只與其他矩陣比較,與自身矩陣之元素比較亦可,例如:

>>A(A<10)=10
A =
100 10 10 100
10  100 100 100
100 100 10 100
10 100 100 10

3.2邏輯操作元

邏輯操作元主要在連繫兩個邏輯值,應用上分為元素間、位元間等兩類。位元間讀者可自行參閱深入書籍,此處稍加介紹元素間之應用。其相關操作元如表3.2。



表3.2 邏輯操作元
操作元範例說明
&C=A&BAB兩元素均為1時為1,否則為0
C=A|B任一A或B為1為正,否則為0
~C=~A非A,即A為0時為1,否則為0
xorC=xor(A,B)僅其中一個為非零時為1,否則為0

       

邏輯操作元常與關係條件諸元配合,作為判斷及分類之標準。下面舉例說明:

>>A=[1 0 1 1 0 0 1], B=[0 0 1 0 1 0 1]
A =
1 0 1 1 0 0 1
B =
0 0 1 0 1 0 1

>>C=A&B
C =
0 0 1 0 0 0 1

>>C=A¦B
C =
1 0 1 1 1 0 1

>>C=~A
C =
0 1 0 0 1 1 0

>>C=xor(A, B)
C =
1 0 0 1 1 0 0

應用時,兩矩陣必須同樣大小,除非其中一個為常數。若元素屬於非零元素(如10,8 -5等等)均視為邏輯1,或為真值。

在實際運作上,前面係直接用&~等三個符號進行邏輯運算,有些地方亦可使用函數型式,即and(A,B)代表A & B;or(A,B)代表A B;not(A)代表 ~A。

利用矩陣本身,進行分類的方式到處可用。這些指令函數亦常使用,其型式如下表:













操作元範例說明
anyC=any(A,dim)A矩陣中之任何行向元素不為0時為1,否則為0。dim=2時為列向。
allC=all(A)矩陣中之任何行向元素均為1時為1,否則為0。dim=2時為列向。
isprimeC=isprime(A)A矩陣中,元素為質素時為1時為1,否則為0
isequalC= isequal(A)僅其中一個為非零時為1,否則為0
isemptyC= isempty(A)若為空矩陣為1,否則為0
isintegerC= isinteger(A)若為整數為1,否則為0
isnumericC= isnumeric(A)若為數值為1,否則為0
isrealC= isreal(A)若為實數為1,否則為0
isfiniteC= isfinite(A)若為定值為1,否則為0
logicalC= islogical(A)將矩陣A轉為邏輯值,非零為1
islogicalC= islogical(A)若為邏輯值為1,否則為0
ischarC= ischar(A)若為文字值為1,否則為0
iscellC= iscell(A)若為細胞陣列為1,否則為0
findC=find(A);C=find(A>B)找尋矩陣A中不為零之序號,存於C找尋矩陣A大於B之序號存於C


例如從前面的原A矩陣找出質因數:

>>A=magic(3)
B=A;
B(~isprime(A))=0

>>A =
8 1 6
3 5 7
4 9 2

B =
0 0 0
3 5 7
0 0 2

在這裡,isprime是尋找矩陣中元素是否為質數,而~isprime即為非質數,將其對應元素設為零即可得到質數之矩陣。這些利用is*之型式逗起來的函數也是屬於矩陣邏輯操作函數,比較常用的有isempty、isequal等等。有興趣的讀者可以參考其輔助檔。

logical則是一個將矩陣轉換為邏輯常數函數,只要零值換為邏輯0值,非零(包括負值)均視為邏輯1。例如:

>> A=-5:5
A =

-5 -4 -3 -2 -1 0 1 2 3 4 5
>> islogical(A)
ans =
0

顯然矩陣A並非邏輯矩陣,但經過如下轉換,即可轉為邏輯矩陣:

>> logical(A)
Warning: Values other than 0 or 1 converted to logical 1.
ans =
1 1 1 1 1 0 1 1 1 1 1

此外,any與 all兩函數則是檢驗一個矩陣是否有非零元素或全為非零元素。兩者均可針對單行或單列檢驗,如:

>>any(A==B)
ans =
1 1 1

顯然它是行向作業,如改為列向則後面加一參數dim=2

>>any(A==B,2)
ans =
0
1
1

注意上述之關係等號必須使用兩個"="號。讀者可自行印證這兩個的結果。在行向上,A與B是沒有完全相等之元素,所以下面用all指令的結果應均為零。

>>all(A==B)
ans =
0 0 0

但由於第二列均為[3 5 7],故下面的結果第二項應為1,其餘均為0。

>>all(A==B, 2)
ans =
0
1
0

有關關係操作元之運用在未來程式寫作上甚為重要。簡潔的關係式可以使程式看起來更為流暢,其運作的時間也會降低。最後find之函數則是尋找某矩陣不為零之元素或具某條件之原始序號,如:

A=fix(rand(3)*10)
A =
9 8 8
9 0 0
4 3 1

>>find(A)'
ans =
1 2 3 4 6 7 9

由結果可知,第五項與第八項為零。此外:

>>find(A>5)'
ans =
1 2 4 7

其結果是:只有第一、二、四、七項元素大於五。若為矩陣間比較,其大小應相同。

3.3 IF條件敘述

條件敘述是一種重要的功能,可以擴大程式判斷的能力。MATLAB最簡單之條件敘述有下面之型式:


if {條件敘述1}
{指令敘述1};
end

if {條件敘述2}
{指令敘述2};
elseif{條件敘述3}
{指令敘述3};
else
{指令敘述4};
end


有關這兩種條件敘述之例子,前面已經有引述。上述條件敘述部份可能為關係式,可能為邏輯判斷值,亦可為邏輯函數的結果,如isreal(A)、isprime(A)等。簡單的邏輯式可以利用關係子給合成為複雜的關係式,但其最後之結果值僅有兩種狀態,即true或 false,或為 0 與1。若結果值不同於1,也視同1考慮。if若成巢狀結構,則必須以end結束。

複合式陳述時,MATLAB並不一定處理全部條件。通常採用部份評估的方式。以if (A&B)為例,只要A的條件評估結果為零時即不再評估B;同理,若為if(AB)時,則只要A為非零之結果,即不再評估B的結果。由於有是項特徵,有些條件陳述可以利用這項功能,如:

if (b~=0) & (a/b>20), …

在此只要(b~=0) 之結果為偽(即b=0),即不會進一步評詁 (a/b>20)這一項,以免造成divide by zero的錯誤。

同理,若條件陳述為

if exist('draw_number.m') & (draw_number(100)>=X)

則第一項為偽時(即檔案draw_number.m不存在時)即不再評估並執行第二項。下面一個例子也是如此,讀者可以自己體會:

if iscell(A) & any(cellfun('isnumeric',A))

下面使用一個抽號的方式作為例子:



function [C]=draw_number(no_of_draw)
% draw ball numbers within ndraw times
%
C=zeros(1,5);
n=1;
while n<=no_of_draw
ball=fix(rand*10);
if ball<2
C(1)=C(1)+1;
elseif ball<4,
C(2)=C(2)+1;
elseif ball<6,
C(3)=C(3)+1;
elseif ball<8,
C(4)=C(4)+1;
else
C(5)=C(5)+1;
end
n=n+1;
end



這是以隨機函數抽號碼球由0至9號分別置於C(1)-C(5)內。看其抽中的號碼分別置入五個箱子中。這個程式完全利用if {}; elseif {};else{};end之型式。最外圍則掛上while…end之迴圈,將在後面討論。執行本程式初期以1000次、5000次及10000次抽到之結果為:

>>[C]=draw_number(1000)
C =
168 209 214 191 218

>>[C]=draw_number(5000)
C =
1010 977 1020 1009 984

>>[C]=draw_number(10000)
C =
2036 2007 1875 2045 2037

依據上項關係得知,抽取到的或然率大概差不多哩!而且即使同樣的抽取次數,其結果亦會有不同,讀者不妨試試。

若以兩矩陣作大小比較時,通常均以元素間之值作比較,但矩陣之大小需相同。下面就實際兩個例子作比較:

>>A=[0 1;2 5], B=[1 0;2 3]
A =
0 1
2 5

B =
1 0
2 3

>>A<B
ans =
1 0
0 0

>>A<B+1
ans =
1 0
1 0

>>A&B
ans =
0 0
1 1

>>A<5
ans =

1 1
1 0

IF之判斷指令可以應用於一般輸入之參數之確定(如前面之evalsum.m程式)或一般變數輸入之交談用。

在此特別強調的一點是,有關"=="與"isequal"這兩個指令之比較。在一般常數間之比較時,此兩個指令之結果應該相同,但在矩陣間比較時,兩個結果可能不同。

>> A==B
ans =
0 0
1 0

>> isequal(A,B)
ans =
0


這種情況也可能發生在其餘類似的指令如 isempty、all、 any等。

3.4 FOR 迴圈之應用

通常for的迴圈都用在重覆執行固定次數的場合,其運算過程中雖亦可利用條件跳出迴圈,但仍以其設定之次數為上限,迴圈每循環一次稱為一回合。其中之expression 通常以純量表示,其型式如k=1:5、K=2:4:20等,等號右邊之最左值為初值,最後值為末值,有中間值時則為每回合之增量。在 MATLAB中 for 迴圈之型式如下:


for var = expression
{statements};
end


例如:



x=0, %先將x矩陣清為零,再進行迴圈
for i = 1:5,
  x(i) = i^2,
end;


x = 0
x = 1
x = 1 4
x = 1 4 9
x = 1 4 9 16
x = 1 4 9 16 25

由這個迴圈之運作,至少可以看出x列向量之大小逐漸因執行過程而擴大的情形。上述迴圈亦可寫成:

>>for i = 1:5, x(i) = i^2, end

但要先記得清除變數x,以免變數有殘餘值。即

>>clear x

若改成遞減的方式亦有同樣的結果:

>>x=0;for i = 5:-1:1, x(i) = i^2, end

x = 0 0 0 0 25
x = 0 0 0 16 25
x = 0 0 9 16 25
x = 0 4 9 16 25
x = 1 4 9 16 25


只是此時之x已先設定為1x5之向量,因為開始時x=5。

使用這種遞圈的觀念亦很容易產生向量矩陣,例如:

>>x=1:1:10
x =
1 2 3 4 5 6 7 8 9 10

這表示產生一個序列向量,起始值為1,增量為1,直到10為止。若增量為1個單位時,上述x向量亦可用x=1:10產生。而增量改為負值時,其順序應相反:

>>x=10:-1:1
x =
10 9 8 7 6 5 4 3 2 1

巢狀迴圈亦常使用,下面是產生修伯特矩陣(Hilbert matrix)的方式:



%hillbert.m for demonstration
%
h=5;
A=zeros(h,h); %預留空間
for k=1:h
for m=1:h
A(k,m)=1/(k + m -1);
end
end
A


>>hillbert
A =
1.0000 0.5000 0.3333 0.2500 0.2000
0.5000 0.3333 0.2500 0.2000 0.1667
0.3333 0.2500 0.2000 0.1667 0.1429
0.2500 0.2000 0.1667 0.1429 0.1250
0.2000 0.1667 0.1429 0.1250 0.1111

實際上,初值亦可為向量,然後依序計算:

>> for i=[2 18 3 1], x=i^3,end

x = 8
x = 5832
x = 27
x = 1


比較深入之運用,for後面之陳述也可以使用矩陣,只是每次以該矩陣之行向為單位,一行一行地執行,如


>>eye(3)

ans =

1 0 0
0 1 0
0 0 1

>> for m=eye(3)*3, m, end
m =
3
0
0
m =
0
3
0
m =
0
0
3


由於eye(3)是一個單位矩陣(對角線元素均為1其餘為0),故每回合所得之m值為其行向量值。在特別情況下,讀者也可以利用這項功能。這與下面之執行結果略有不同。


>> m=eye(3)*3
m =
3 0 0
0 3 0
0 0 3

3.5 While 迴圈

while與for最大不同點在於前者係依某種條件結束迴圈,故無法如for執行時可預知其迴圈數目。While 迴圈之其型式如下:



while expression
{statements};
end




在while後面之expression只要其評估結果為真,此迴圈會繼續執行,因此無法確知其執行次數,為此有時必須另設計數器計算次數。有時為使其形成永久迴圈,亦可採用while 1,…end之型式,使迴圈永遠有效。此時必須有特殊狀況才能跳出,或停電關機才能終止。前面之例子中,已有示範如何使用while...end這個迴圈。茲再以計算前N個自然數和之例說明:


% cumsumx.m for demonstration
N = 50;
n = 0;
sum = 0; %先設定儲存變數
while N>=n, %只要n 小於50,下面之指令陳述繼續執行
n = n+1;
sum = sum+n;
end
sum

sum=

1275


上面程式執行後,其值應為1275。此小程式中可以得知大寫的N與小寫的n是兩個不相干的變數。故大小寫的變數不能混淆,否則會有很大的麻煩。

下面的範例為利用while…end計算一台機器所能處理之最小值。這個最小值因電腦軟體不同而異,在MATLAB中有一個常數eps是直接可以拿來使用,如:

>>eps

eps=

2.2204e-016

此處用程式亦可計算該值:


% epsilon.m for demonstration
eps =1;
while (1+eps) >1
eps=eps/2;
end
eps=eps*2;

>>epsilon
eps = 2.2204e-016

其結果相同。此處之while…end則是評估(1+eps) >1之真偽決定結果,直到eps小到電腦無法分辨其大小為止。

3.6 SWITCH條件式

與前面之if-else-end 敘述相同的條件式,有時亦可利用switch 之敘述方式,使結構更條理化。switch 之陳述型式如下:

  
switch switch_expr
case case_expr,
statement, ..., statement
case {case_expr1, case_expr2, case_expr3,...}
statement, ..., statement
...
otherwise,
statement, ..., statement
end


在這裡只要第一個 CASE 的條件相符合時,其後之敘述指令即會被執行。若CASE 之條件是以數個單元之集合表示(如第二個CASE),則只要其中一單元符合SWITCH後之條件,即會執行第二CASE下之指令陳述。若所有條件均不符合,則會執行OTHERWISE後之指令陳述。不論如何,僅其中一個CASE的內容會被執行,執行後即會跳到 END後面之敘述去執行下一段之指令。這種方式可以說是最結構化之條件式,也最常用。

在沒離開這個題目之前,仍有幾句話附帶說明,在SWITCH後面之敘述switch_expr,可以為常數或字串。故只要switch_expr==case_expr 為真,則該CASE下之陳述會被執行。

下面為利用前面if的例子draw_number程式改寫為switch case之應用:

function [C]=draw_number2(no_of_draw)
% draw ball numbers within ndraw times
% Using switch case
C=zeros(1,5);
n=1;
while n<=no_of_draw
ball=fix(rand*10);
switch ball
case {0,1}
C(1)=C(1)+1;
case {2,3}
C(2)=C(2)+1;
case {4,5}
C(3)=C(3)+1;
case {6,7}
C(4)=C(4)+1;
otherwise
C(5)=C(5)+1;
end
n=n+1;
end


執行結果:

>>draw_number2(1000)
ans =
186 197 211 199 207

>>draw_number2(5000)
ans =
1017 1028 934 999 1022

>>draw_number2(10000)

ans =
2014 2021 1908 2036 2021

執行結果大致相同。switch case比較適合列舉型選項,if指令則可適用於各種範圍的選擇,其應用均很廣泛。

3.7 BREAK指令

break這一個指令常與迴圈while 或for 配合使用,可以自迴圈中跳出至上一層之迴圈,在巢狀迴圈中,break僅能自最內圈跳出圈外,到其對應之end後繼續執行。下面之範例配合while 1…end執行輸入的問答,只有回答N才能利用break跳出迴圈。


% Using break to terminate the execttion.
n=1;
while 1
a=input('Continue? (Y/N) ','s');
A=upper(a);
if isempty(A), A='Y';end
if A=='N', break; end
if A=='Y',
disp('OK, you know how to use BREAK, don''t you?')
disp(['You just try ',int2str(n),' time(s)! Do it again.'])
n=n+1;
else
disp('I don''t understand, please input again.')
end
end


執行之結果:

>>show_break
Continue? (Y/N) d
I don't understand, please input again.
Continue? (Y/N) y
OK, you know how to use BREAK, don't you?
You just try 1 time(s)! Do it again.
Continue? (Y/N) y
OK, you know how to use BREAK, don't you?
You just try 2 time(s)! Do it again.
Continue? (Y/N) n

3.8 CONTINUE & RETURN

continue與break兩指令具有相同的功能,但break大部份都是跳到外圈,continue則只跳到對應之end後繼續執行,在巢狀if…end內,其功能大致相同。

將上例稍作修正,可以相互比較其功用之不同。本例仍然有使用break中斷執行while…end。所以兩者有不同的功能。


% Using continue to terminate the execttion.
n=1;
while 1
a=input('To proceed? (Y/N) ','s');
A=upper(a);
if isempty(A), A='Y';end
if A=='N', break; end
if A=='Y',
disp('OK, you know how to use CONTINUE, don''t you?')
disp('This statement uses CONTINUE to bypass the rest ones.')
disp(['You have tried ',int2str(n),' time(s)! Try it again.'])
n=n+1;
continue
end
disp('I don''t understand, please input again.')
end


>> show_continue
To proceed? (Y/N) y
OK, you know how to use CONTINUE, don't you?
This statement uses CONTINUE to bypass the rest ones.
You have tried 1 time(s)! Try it again.
To proceed? (Y/N) n

return也是中斷指令,但通常用於被呼叫之函數中,一旦遇到return指令,執行權將回到呼叫程式繼續執行。


function dd=show_det(A)
% Find the determinant of A
if isempty(A),
dd=1;
return
else
dd=det(A);
end


執行例:

>> d=show_det('')
d = 1
>> d=show_det([1 2 3;4 -5 6;7 2 1])
d = 188

3.9 模組化設計

一般而言,程式均可以配合流程控制結構撰寫出所需的程式內容,而Matlab與其他電腦一樣,具備有結構化程式設計之能力。通常在這些程式設計中,不再使用goto這個指令,雖然這在Fortran與 Basic語言中,是常用的指令,但這種指令碼很容易產生秩序混亂的程式設計,容易使流程糾結不清。

一般結構化程式具有下列特色:

  1. 容易撰寫:先就問題內容瞭解後,再處理細節,較容易產生簡潔的大綱。
  2. 可重複使用模組:具有相同功能之程式模組可以一再使用。
  3. 容易除錯:流程簡單,容易除錯。
  4. 可以合作開發:在研發過程中可以將整個程式進行切割,最後再進行整合。
  5. 容易辨識及瞭解:程式寫作完成後,有如一篇完整脈絡的報告,容易按圖索驥。

結構化程式設計通常採用由上而下的設計,開始時才提綱要,然後再逐步進入細節。在文件說明上更容易利用結構圖及流程圖表示。前者為描述程式部份如何相互聯結,因此可以很清楚地顯示整個程式之架構。

例如一套乾燥系統之程式可以將整個程式作如下之分類:

  主程式
     。濕氣特性程式
。乾燥理論程式
     。乾燥機設備程式
。進料程式
。循環程式
。出料程式


流程圖對於特定程式之執行、判斷及分支之狀況有比較清晰的描述。如此,任何對程式有興趣的人,就流程之變換可以立即瞭解轉折點的位置及程式進行時所採取的步驟。結構圖及流程圖對解題上有也相當大的助益,在偵錯過程中,更能一目了然。





















在程式撰寫過程中,將程式文作作適當說明也是一件相當重要的步驟。對後來要對程式有深入瞭解時有相當大的幫助。在MATLAB中,可以充分利用%這個控制碼,因為在控制碼%後的文字程式並不執行,故可以用這個方式進行說明。在指令後面亦可加%進行個別說明,作為程式之註解。在程式文件中,亦可考慮下列動作,以協助程式清晰度:

  1. 為程式之函數、變數選取適當名稱,以反應該變數的量與意義。
  2. 在程式當中多加註解,並作變數之意義詳加定義。
  3. 配合結構圖加以說明及分區歸類。
  4. 配合流程圖說明判別功能點。
  5. 使用虛擬碼(Pseudocode)協助釐清程式的細節。配合文字及演算公式,使計算步驟清楚顯現。


有關於這些細節之運用將在未來章節中配合說明。