星期四, 12月 28, 2006

Topic: 人體受外力牽拉後引發之姿勢反應_肌肉活化情形_rev.2

程式更新描述:
1. 加入各個牽拉強度下3次測試IEMG平均值及標準差計算,並將每位受試者各個測試以及各個牽拉強度(平均值與標準差)的IEMG計算結果輸出於同一檔案中,並以.xls檔案形式儲存於特定資料夾中。
2. Missing data的處理
(1)當原始EMG資料中任一塊肌肉訊號有所遺漏時,其計算出來的IEMG數值會以NaN顯示,程式可判斷'如果同一次測試中,有任一塊肌肉IEMG數值缺乏時,則該次測試所有肌肉的IEMG數值以-99表示,且皆不列入最後的平均值與標準差之計算'。
(2)另外,程式亦可判斷'當某次測試的原始EMG資料檔案不存在時,則該次測試所有肌肉的IEMG數值以-99表示,且皆不列入最後的平均值與標準差之計算'。

程式內容:

% IEMG6.m
% function: IEMG calculation (including mean and std of the value of each
% perturbation amplitude) for all trials of each subject
% onset data were retrived automatically
% if the file of raw data is not existed or the data of the measurments of a trial is not recorded, print all the IEMG data of that trial '-99'
% the IEMG results were output and saved as .xls
% raw_data: 1_Time(ms) 2_Fy 3_Mx 4_MGr 5_HAMr 6_PARAr 7_GMr 8_Pullonset
% range from (600ms prior to the pull onset) to (5000ms after the pull
% onset)
% baseline EMG:time interval between (500ms prior to the pull onset) to (150ms priot to the pull onset)
% IEMG:time interval between (150ms after the burst onset) to (500ms after the burst onset)

global EMGdata_path;
EMGdata_path = input('Please input the folder of the EMG data: [e.g.: D:\\Physical Therapy\\Graduate School\\Thesis\\Pilot study for Thesis\\Data\\Processed data for Matlab]','s');

if isempty(EMGdata_path)
EMGdata_path = 'D:\Physical Therapy\Graduate School\Thesis\Pilot study for Thesis\Data\Processed data for Matlab';
end

cd(EMGdata_path);
Subject_information = textread('Subject_information.txt','%s');
Subject_amount = size(Subject_information,1);

EMGoutput_path = ['D:\Physical Therapy\Graduate School\Thesis\Pilot study for Thesis\Results\Data output from Matlab'];
if ~exist(EMGoutput_path,'dir')
mkdir('D:\Physical Therapy\Graduate School\Thesis\Pilot study for Thesis\Results\Data output from Matlab'); % creates the directory dirname in the current directory, if dirname represents a relative path.
end

for i = 1:Subject_amount

subjectoutput_folder = ['D:\Physical Therapy\Graduate School\Thesis\Pilot study for Thesis\Results\Data output from Matlab','\',Subject_information{i}];
if ~exist(subjectoutput_folder,'dir') % dir 顯示目前目錄下的檔案名稱
mkdir('D:\Physical Therapy\Graduate School\Thesis\Pilot study for Thesis\Results\Data output from Matlab',Subject_information{i});
end

Onsetdata_xls = xlsread([EMGdata_path,'\',Subject_information{i}(1:4),'\','Time events for EMG and pullonset signials.xls'],'B3:F14'); % [pullonset,MGronset,HAMronset,PARAronset,GMronset]

data_name=0;
Onsetdata=0;
IEMG=0;

fprintf('\n %s\t\t rMG (V*ms)\t rHAM (V*ms)\t rPARA (V*ms)\t rGM (V*ms)\n',Subject_information{i}(1:4))

IEMG_fid = fopen([subjectoutput_folder,'\',Subject_information{i},'Integrated EMG','.xls'],'w');
fprintf(IEMG_fid,' \t rMG (V*ms)\t rHAM (V*ms)\t rPARA (V*ms)\t rGM (V*ms)\n');

for j = 1:4 % 4 perturbation amplitudes
for k = 1:3 % 3 trials
data_name = [Subject_information{i}(1:4),'C',num2str(1),'V',num2str(j),'R',num2str(k)];
Onsetdata = Onsetdata_xls((j-1)*3+k,:);
rawdata_file = [EMGdata_path,'\',Subject_information{i}(1:4),'\',data_name,'.xls'];
if exist(rawdata_file)
raw_data = xlsread([EMGdata_path,'\',Subject_information{i}(1:4),'\',data_name,'.xls'],'C19:J5618');
baseline = mean(raw_data(101:451,4:7));

rawIEMG(1) = (1/2*(raw_data((Onsetdata(2)-Onsetdata(1)+149),4)+raw_data((Onsetdata(2)-Onsetdata(1)+499),4))+sum(raw_data((Onsetdata(2)-Onsetdata(1)+150):(Onsetdata(2)-Onsetdata(1)+498),4))); % unit:V*ms
rawIEMG(2) = (1/2*(raw_data((Onsetdata(3)-Onsetdata(1)+149),5)+raw_data((Onsetdata(3)-Onsetdata(1)+499),5))+sum(raw_data((Onsetdata(3)-Onsetdata(1)+150):(Onsetdata(3)-Onsetdata(1)+498),5))); % unit:V*ms
rawIEMG(3) = (1/2*(raw_data((Onsetdata(4)-Onsetdata(1)+149),6)+raw_data((Onsetdata(4)-Onsetdata(1)+499),6))+sum(raw_data((Onsetdata(4)-Onsetdata(1)+150):(Onsetdata(4)-Onsetdata(1)+498),6))); % unit:V*ms
rawIEMG(4) = (1/2*(raw_data((Onsetdata(5)-Onsetdata(1)+149),7)+raw_data((Onsetdata(5)-Onsetdata(1)+499),7))+sum(raw_data((Onsetdata(5)-Onsetdata(1)+150):(Onsetdata(5)-Onsetdata(1)+498),7))); % unit:V*ms

IEMG = rawIEMG-baseline; % [rMG rHAM rPARA rGM]
if ~isempty(find(isnan(IEMG))) % 如果有其中一個參數出現missing,則那一個trial內所有的參數皆不列入最後的平均計算
IEMG(1,1:4) = -99;
end
else
IEMG(1,1:4) = -99;
end

IEMG_command = ['IEMG_cell(',num2str(k),',:) ={IEMG};'];
eval(IEMG_command)

fprintf(' %s\t %f\t %f\t %f\t %f\n',data_name,IEMG)

fprintf(IEMG_fid,' %s\t %.3f\t %.3f\t %.3f\t %.3f\n',data_name,IEMG);
end % for k loop

IEMG_array = reshape([IEMG_cell{:}],4,3)';
for l = 1:4
if ~isempty((find(IEMG_array(:,l) ~= -99)))
mean_IEMG(1,l) = mean(IEMG_array(find(IEMG_array(:,l) ~= -99),l));
std_IEMG(1,l) = std(IEMG_array(find(IEMG_array(:,l) ~= -99),l));
else
mean_IEMG(1,1:4) = -99;
std_IEMG(1,1:4) = -99;
end
end

fprintf(IEMG_fid,' \n');
fprintf(IEMG_fid,' %s\t %.3f\t %.3f\t %.3f\t %.3f\n','Mean',mean_IEMG);
fprintf(IEMG_fid,' %s\t %.3f\t %.3f\t %.3f\t %.3f\n','Std',std_IEMG);
fprintf(IEMG_fid,' \n');
end % for j loop

fclose(IEMG_fid);

end % for i loop

討論:
待完成事項:
1. 將程式一開始改為3種mode選擇:mode1-計算IEMG;mode2-單純繪出代表的EMG圖形;mode3-計算肌肉活化程度(IEMG)與干擾強度之相關性並繪圖。

星期六, 12月 23, 2006

Blog 小程式相關網站

@ 計時器

※ ClockLink
http://www.clocklink.com/
有各種不同分類的計時器,不過相較於日本人做的版本就有點死板的感覺

※ブログパーツの森
http://www.blogpartsnomori.com/
提供不同款式的小熊造型時鐘,超可愛 ><

※ポストペットブログアクセサリ(PostPet Blog Accessory) http://www.postpet.so-net.ne.jp/webmail/blog/ 多款造型可愛的寵物小時鐘,現在網頁上的企鵝計時器就是從這裡找到的喔....

@ 另類小程式

※ Moon Phases
http://www.calculatorcat.com/moon_phases/moon_phases.phtml
能夠知道現在月亮的陰晴圓缺喔....

@ 提供多種Blog小程式

※ Blog Apartments
http://blog-apart.com/index.html
分為三大類:Utility Blog Parts、Game Blog Parts,還有Interactive Blog Parts

※ブログパーツ.com
http://www.blog-parts.com/
分成十六個項目,包括計時器、溫度計、占卜、天氣預報等小程式。有些需要先登錄才能下載喔。

星期四, 12月 21, 2006

matlab課堂筆記 061222_第十一章 統計與迴歸

相關係數 corrcoef

R = corrcoef(X)
R = corrcoef(x,y)
[R,P]=corrcoef(...)


R代表correlation coefficient;P即代表判斷是否有統計意義的p

例題:
>> x=1:5;y=x.^3;
>> [R,P]=corrcoef(x,y)

R =
1.0000 0.9431
0.9431 1.0000
P =
1.0000 0.0161
0.0161 1.0000
>> [i,j] = find(P<0.05)

i =
2
1
j =
1
2

R(1,1)代表x和x之間的相關係數 → 一定為1
R(2,1)代表y和x之間的相關係數 → 一定為1
R(1,2)代表x和y之間的相關係數
R(2,2)代表y和y之間的相關係數


群組函數 grpstats

應用於計算複雜分組的平均問題。

means = grpstats(X, group)
[means, sem, counts, name] = grpstats(X, group, whichstats)
grpstats(x, group, alpha)


輸入參數中X為求平均值之對象(可為多行,其平均結果也會多行)。group則為與X同列長之陣列,可能由多項分組之向量組成,其內容可為字串列或細胞陣列之文字,如{G1 G2 G3}。在輸出項中,第一項means為群組平均,sem為組內標準差,counts為各組之項數,name則為各組之名稱。

例題:
x =
1 2 3 4 5 6 7 8 9
>> group=[1 1 1 1 2 2 2 2 2]
group =
1 1 1 1 2 2 2 2 2
>> [m,s,c]=grpstats(x,group)
m =
2.5000
7.0000
s =
0.6455
0.7071
c =
4
5

特殊指令 normrnd

normrnd(m,s)

可由一個平均值為m,標準差為s的平均分佈中隨機產生一個數值。

堆疊矩陣 repmat

B = repmat(A,m,n)
B = repmat(A,[m n])
B = repmat(A,[m n p...])


以A之內容堆疊在一(M x N)的矩陣B中。B矩陣之大小由MXN及A矩陣之內容決定。

例題:可以repmat指令作出與ones(3,3)*5相同的結果
>> repmat(5,3,3)

ans =

5 5 5
5 5 5
5 5 5

和使用ones(m,m)*n的方式相比的優點是
使用repmat指令可以在A置放文字串!!!

Topic: 人體受外力牽拉後引發之姿勢反應_肌肉活化情形_rev.1

程式更新描述:
1. 將onset時間輸入方式改為由程式直接從已建立之excel檔案內擷取所需資料,不須經由對話窗輸入。
2. 新程式可同時讀取不同資料夾(依受試者ID不同做區別)內各個測試原始EMG data,並分別與該測試的onset時間資料相對應,如此可同時計算多位受試者以及多個測試的IEMG結果。
3. 將計算結果依每位受試者的測試順序做排序,俾輸出於指令窗中。

注意事項:
1. 實驗結果可儲存於任一路徑,若無指定路徑,則依預設路徑(D:\Physical Therapy\Graduate School\Thesis\Pilot study for Thesis\Data\Processed data for Matlab)進行原始資料之讀取。預設路徑內除了原始的實驗結果之外,亦須包括一Subject_information.txt檔案,紀錄欲分析之受試者ID,若欲分析不同受試者,則須由此檔案做變更。
2. 各個受試者的實驗結果以.xls檔案形式儲存於個別資料夾中(命名方式如,He01、He02、St05等;He/St代表不同受試者組別),各個受試者之資料夾中須包括所有測試的原始EMG data(命名方式如,He01C1V1R1、He02C1V2R3等;V代表不同牽拉干擾強度,R代表不同測試)以及onset時間檔案(Time events for EMG and pullonset signials.xls)

程式內容:

% IEMG4.m
% function: IEMG calculation for all trials of each subject
% onset data were retrived automatically
% the IEMG results were displayed on the command window
% raw_data: 1_Time(ms) 2_Fy 3_Mx 4_MGr 5_HAMr 6_PARAr 7_GMr 8_Pullonset
% range from (600ms prior to the pull onset) to (5000ms after the pull onset)
% baseline EMG:time interval between (500ms prior to the pull onset) to (150ms priot to the pull onset)
% IEMG:time interval between (150ms after the burst onset) to (500ms after the burst onset)

global EMGdata_path;
EMGdata_path = input('Please input the folder of the EMG data: [e.g.: D:\\Physical Therapy\\Graduate School\\Thesis\\Pilot study for Thesis\\Data\\Processed data for Matlab]','s');

if isempty(EMGdata_path)
EMGdata_path = 'D:\Physical Therapy\Graduate School\Thesis\Pilot study for Thesis\Data\Processed data for Matlab';
end

cd(EMGdata_path);
Subject_information = textread('Subject_information.txt','%s');
Subject_amount = size(Subject_information,1);


for i = 1:Subject_amount


Onsetdata_xls = xlsread([EMGdata_path,'\',Subject_information{i}(1:4),'\','Time events for EMG and pullonset signials.xls'],'B3:F14'); % [pullonset,MGronset,HAMronset,PARAronset,GMronset]
data_name=0;
Onsetdata=0;
IEMG=0;

fprintf('\n %s\n',Subject_information{i}(1:4))

for j = 1:4 % 4 perturbation amplitudes
for k = 1:3 % 3 trials
data_name = [Subject_information{i}(1:4),'C',num2str(1),'V',num2str(j),'R',num2str(k)];
if exist([EMGdata_path,'\',Subject_information{i}(1:4),'\',data_name,'.xls'])
data_name = [Subject_information{i}(1:4),'C',num2str(1),'V',num2str(j),'R',num2str(k)];
else
data_name = [Subject_information{i}(1:4),'C',num2str(1),'V',num2str(j),'R',num2str(k),'r'];
end
Onsetdata = Onsetdata_xls((j-1)*3+k,:);
raw_data = xlsread([EMGdata_path,'\',Subject_information{i}(1:4),'\',data_name,'.xls'],'C19:J5618');
baseline = mean(raw_data(101:451,4:7));

rawIEMG(1) = (1/2*(raw_data((Onsetdata(2)-Onsetdata(1)+149),4)+raw_data((Onsetdata(2)-Onsetdata(1)+499),4))+sum(raw_data((Onsetdata(2)-Onsetdata(1)+150):(Onsetdata(2)-Onsetdata(1)+498),4))); % unit:V*ms
rawIEMG(2) = (1/2*(raw_data((Onsetdata(3)-Onsetdata(1)+149),5)+raw_data((Onsetdata(3)-Onsetdata(1)+499),5))+sum(raw_data((Onsetdata(3)-Onsetdata(1)+150):(Onsetdata(3)-Onsetdata(1)+498),5))); % unit:V*ms
rawIEMG(3) = (1/2*(raw_data((Onsetdata(4)-Onsetdata(1)+149),6)+raw_data((Onsetdata(4)-Onsetdata(1)+499),6))+sum(raw_data((Onsetdata(4)-Onsetdata(1)+150):(Onsetdata(4)-Onsetdata(1)+498),6))); % unit:V*ms
rawIEMG(4) = (1/2*(raw_data((Onsetdata(5)-Onsetdata(1)+149),7)+raw_data((Onsetdata(5)-Onsetdata(1)+499),7))+sum(raw_data((Onsetdata(5)-Onsetdata(1)+150):(Onsetdata(5)-Onsetdata(1)+498),7))); % unit:V*ms

IEMG= rawIEMG-baseline; % [rMG rHAM rPARA rGM]

fprintf(' %s\t %f\t %f\t %f\t %f\n',data_name,IEMG(1),IEMG(2),IEMG(3),IEMG(4))
end % for k loop
end % for j loop
end % for i loop

執行結果:



討論:
1. 目前程式內容,僅限於將計算結果輸出於指令窗中,待完成部分包括:(1)將各個trial之IEMG計算值和每3個trial之平均值和標準差(代表各個牽拉強度)之執行結果以excel檔案形式輸出並儲存於特定資料夾中。
(2) 選出欲代表的EMG活化圖形之trial(以對話方式進行),並將圖形輸出且儲存於既定資料夾內

星期四, 12月 14, 2006

matlab課堂筆記 061215_內插法

內插法 interp1
當查表時,如果欲查的值不在表格內,可用interp1求得
舉例:

>> x=0:10;y=sin(x);[x' y']

ans =
0 0
1.0000 0.8415
2.0000 0.9093
3.0000 0.1411
4.0000 -0.7568
5.0000 -0.9589
6.0000 -0.2794
7.0000 0.6570
8.0000 0.9894
9.0000 0.4121
10.0000 -0.5440

>> y1=interp1(x,y,4.5)

y1 =
-0.8579

也可以由y值反求x值
>> x1=interp1(y,x,0.8)
x1 =
2.3488

Topic: 人體受外力牽拉後引發之姿勢反應_肌肉活化情形



(Horak et al., 1992)






前言:
當人體在站立姿勢下接受一外在的平衡干擾時,下肢與軀幹肌肉的活化程度可以反映人體受到外在干擾後引發之姿勢反應。而當人體在站姿下受到外力干擾而產生前向的身體擺盪時,背側的肌肉主要負責活化以維持或回復平衡,避免跌倒發生。過去文獻指出,當站立於一塊會產生水平後向移動之力板上時,隨著干擾強度增加(力板移動速度加快或移動位移增加),受試者肌肉活化程度亦會隨之增加。但,由腳底板所接受之干擾與從腰間承受之外在干擾對於受試者的姿勢干擾程度與形式是有所不同的,因此,當受試者由腰間接受干擾後,其肌肉活化情形是否亦會隨著干擾強度增加而增加則值得探討。


實驗目的:
探討人體在站姿下從腰部位置受到之前向牽拉干擾時背側肌肉之活化情形,且比較在不同牽拉強度干擾下,人體反應是否會隨之變化。



實驗流程:
在受試者肚臍位置穿戴一腰帶,將三牽拉軸之牽拉線固定於受試者腰帶上(如上圖),由三牽拉軸產生之合力提供為干擾受試者站立平衡之牽拉力量。受試者在受到牽拉之後必須在不跨步的前提下盡可能地維持平衡不跌倒。牽拉強度共有4種 (固定牽拉速度為9cm/s,牽拉位移則依序增加:3,5,7,9cm),每個牽拉強度情境各進行3次測試。以肌電圖所量測到的肌肉活化情形作為受試者受牽拉後的姿勢反應依據,受測肌肉包括:腓腸肌 (MG)、膕旁肌(HAM)、臀中肌(GM),以及脊旁肌(PARA)。



牽拉位移 (cm)牽拉速度 (cm/s)
V139
V259
V379
V499


MATLAB程式目的:
1. 由既訂之檔案途徑擷取資料:一為原始檔案(包括時間與肌肉活化大小(以voltage表示)),一為選取後之onset時間(牽拉開始以及各肌肉開始活化時間)。
2.
功能一(mode 1):由原始EMG資料(excel檔案)計算不同牽拉強度所引發之integral EMG大小(IEMG計算方式為,肌肉開始活化後150~500毫秒減去牽拉開始前150~500毫秒時間內之IEMG大小做標準化)。並計算各個牽拉強度下各受測肌肉IEMG之平均值與標準差,以excel格式輸出並儲存於既定之資料夾內。
功能二(mode 2):畫出各個牽拉強度下EMG活化情形(各以其中一個具代表性的trial呈現),並儲存於既定之資料夾內。
Flowchart


以下為部分完成之程式內容
程式內容:
% IEMG.m

% function:IEMG calculation and EMG figure for one trial of a subject

% onset data were inputted by key-in

% IEMG results and EMG figure were displayed on the command window% raw_data: 1_Time(ms) 2_Fy 3_Mx 4_MGr 5_HAMr 6_PARAr 7_GMr 8_Pullonset

% range from (600ms prior to the pull onset) to (5000ms after the pull onset)

% baseline EMG:time interval between (500ms prior to the pull onset) to (150ms priot to the pull onset)

% IEMG:time interval between (150ms after the burst onset) to (500ms after the burst onset)


pullonset=input('please input the onset time of the perturbation (ms):');

MGronset=input('please input the onset time of the Medial Gastrocnemius muscle activation (ms):');

HAMronset=input('please input the onset time of the Hamstring muscle activation (ms):');

PARAronset=input('please input the onset time of the Paraspinal muscle activation (ms):');

GMronset=input('please input the onset time of the Gluteus Medius muscle activation (ms):');

raw_data=xlsread('D:\Physical Therapy\Graduate School\Thesis\Pilot study for Thesis\Data\ASCII output\Selected channels_Fy_Mx_rMG_rHAM_rPARA_rGM_pull onset\He01C1V1R1.xls','C19:J5618');

EMGonset=[MGronset,HAMronset,PARAronset,GMronset];

baseline=mean(raw_data(101:451,4:7));

for i=1:4

rawIEMG(i)=(1/2*(raw_data((EMGonset(i)-pullonset+149),i+3)+raw_data((EMGonset(i)-pullonset+499),i+3))+sum(raw_data((EMGonset(i)-pullonset+150):(EMGonset(i)-pullonset+498),i+3))); % unit:V*ms

end

IEMG=rawIEMG-baseline;

IEMG_MGr=IEMG(1),IEMG_HAMr=IEMG(2),IEMG_PARAr=IEMG(3),IEMG_GMr=IEMG(4)



% 以下為繪圖指令

muscle_name = {'MG (V)' 'HAM (V)' 'PARA (V)' 'GM (V)'};

linecolor=['b' 'r' 'g' 'm'];

for i=1:4

subplot(4,1,i),fig_EMG=plot((raw_data(:,1)),(raw_data(:,i+3)),linecolor(i));xlabel('Time (ms)');box off;hold on;onsetline=plot(pullonset*ones(100,1),0:max(raw_data(:,i+3))/99:max(raw_data(:,i+3)),'k');text(pullonset,max(raw_data(:,i+3)),'pull onset');

ylabel_command = ['ylabel(''',muscle_name{i},''')'];

eval(ylabel_command)

end

執行結果:

IEMG數值


EMG圖形


討論:
1. 實驗結果的蒐集會依實驗情境不同而讓受測者進行多個測試,所以結果大多會以多個檔案分別儲存(以測試名稱做區別)。程式的目的主要是想將實驗結果以最簡便的方式做運算,不須經過太多人為的介入,所以程式的重點將大部分著重在資料的處理部份,減少對話框方式的進行,讓使用者可以在同一時間內處理多筆資料。
2. 目前的程式內容只單純建立IEMG計算方式以及EMG活化之圖形呈現格式,且牽拉和肌肉活化onset時間是以對話方式呈現。之後待完成的部分還包括:
(1) 將onset時間輸入方式改為由程式直接從已建立之excel檔案內擷取所需資料。
(2) 繼續完成在同時間內讀取多個trial之測試結果,並將各個trial之IEMG計算值和每3個trial之平均值和標準差(代表各個牽拉強度)之執行結果以excel檔案形式輸出並儲存。
(3) 選出欲代表的EMG活化圖形之trial(以對話方式進行),並將圖形輸出且儲存於既定資料夾內

星期五, 12月 08, 2006

matlab課堂筆記061207_稀疏矩陣

稀疏矩陣的應用 sparse精簡一般含有零元素甚多的矩陣,僅就非零元素之位置加以陳述。可以將一個含有許多零的元素矩陣改寫為位址表示,但其實質仍如原來之矩陣,功能相同 (其大小仍然與原來之M矩陣相同,而且可以直接與其他矩陣相互運算)。可使用full指令復原。
S = sparse(A) (A為矩陣之變數)
S = full(A)

稀疏函數 spfun
f = spfun(@fun,S)
對稀疏矩陣S之非零元素作函數(@fun)取代

星期二, 12月 05, 2006

MATLAB課堂作業:學生期中考成績 (敘述統計_資料分佈)

前言:
敘述統計(descriptive statistic)能夠初步呈現整體資料情形,透過直方圖圖形的表示方法讓人更能清楚地了解資料之分佈情況。以下,以"學生期中考成績"做為範例。

老師網頁程式如下:
原程式

程式指令:

% statistic.m
% Degree distribution
% Data Input

n=0;
while 1
A=input('輸入重量資料(按enter輸入結束): '); % '請輸入學生期中考成績(按enter輸入結束):'
if isempty(A), break;end
n=n+1;Y(n)=A;
end
%
disp(['==>總計輸入 ' num2str(n) ' 筆資料.'])
disp(['最小值:' num2str(min(Y)) ' 最大值:' num2str(max(Y))])
fprintf('平均值:%0.2f 標準差:%0.2f\n',mean(Y),std(Y))
% 加入平均值和標準差的計算,並且設定數值的格式為小數點後兩位,所以改用fprintf指令
while 1
Rk=input('輸入初始點及級矩(以空格分開)','s');
if ~isempty(Rk),
rank=str2num(Rk);
r1=rank(1);r2=rank(2);
break; % 如果按沒有輸入直接enter的畫會重複出現問句
end
end
D=zeros(fix((max(Y)-min(Y))/r2)+2,1); % fix指令為捨去小數部分,D矩陣為結果中的frequency部分。原本只須加1,但是因為結果中有<的部分所以加2
for i=1:n
if Y(i)<r1,
D(1)=D(1)+1;
else
J=fix((Y(i)-r1)/r2)+2;D(J)=D(J)+1;
end
end
disp([' <' num2str(r1) ' ' num2str(D(1))]);
for i=2:length(D)
disp([num2str(r1+(i-2)*r2) '-'...
num2str(r1+((i-1)*r2)) ' ' num2str(D(i))])
end
% 以bar指令畫出資料分佈頻率之圖形 (長條圖)
H=zeros(fix((max(Y)-min(Y))/r2)+2,1);
H(1)=r1-(r2)/2;
for j=2:length(D)
H(j)=r1+(j-2)*r2+(r2)/2;
end
bar(H,D)
xlabel('Value') % '期中考分數'
ylabel('Frequency') % '學生個數'
title('Distribution of the Weight')% 'T大xx系學生期中考成績分佈圖'

執行結果:
統計結果
應用:學生期中考成績

討論:
1. 加入平均值和標準差資料更能清楚呈現資料。
2. 此程式中,若輸入值皆為正值,而起始點輸入為0時,程式執行會出現錯誤,這一部分需再修改。

星期五, 12月 01, 2006

matlab課堂筆記 061201_第十章 聯立方程式解

行列式與反矩陣
AA-1=I
A-1A=I
I稱為單位矩陣(Identity matrix),為一方矩陣,且其對角線元素值均為1。以eye指令可製作此矩陣。
矩陣A之反矩陣(A-1)以inv(A)表示。A矩陣之行列式值需不得為零,否則其反矩陣不能存在 (以det(A)作測試,計算其行列式值,若det(A)=0,則矩陣A為奇異矩陣)。

線性聯立方程式矩陣解法
聯立線性方程式可以化成AX=C之型式:
3x +4y =10
5x -2y =8

(1) 左除法求解:利用MATLAB求解時,只要用矩陣倒除即可,或稱為左除法。即X=A\C
A=[3 4;5 -2];
C=[10;8];
X=A\C
X =
2 1
(2) 反矩陣法求解:A-1AX=A-1C  或A-1AX=IX=X=A-1C
A=[3 4;5 -2];
C=[10;8];
X=inv(A)*C
X =
2.0000 1.0000
答案和用左除法所求得的解相同,但是左除法較常見。

星期四, 11月 30, 2006

matlab課堂筆記 061201_第九章 繪圖 II

直方圖 hist
大多應用於統計學上處理大量資料之分佈頻率。
N=hist(Y,m)
將Y資料之分佈範圍區分為m等分組,然後依資料之大小決定其落點,注意此時m僅為常數。若沒有輸入m項,則其預定組數為10。


繪圖指令 xtick xticklabel

例題:

set(gca,'xtick',1:n,'xticklabel',year)

gca為宣告前一個畫出來的圖形
xtick為設定x軸的間距
xticklabel為設定x軸label的名稱

星期三, 11月 29, 2006

MATLAB課堂作業: 作業九

1. What's difference between the commands "line" and "plot", explain by using your own examples.
ANS:
line指令與plot指令的差別主要如下:
(1) 以line指令可在同一個圖執行多次線段之繪製,以多行對應表示,每一行代表一線;但是使用plot指令時,於圖形視窗內新產生的圖形會覆蓋掉原有的圖形,如果想在已存在的圖中加入新圖,則需要在繪圖指令中間加上hold on指令,使新產生的圖形疊加在原有圖形上面。

line 指令
>> x=1:1000;
>> y1=x.^(1/2);y2=x.^(1/3);y3=x.^(1/4);
>> h1=line(x,y1);
>> h2=line(x,y2);
>> h3=line(x,y3);

plot 指令
>> x=1:1000;
>> y1=x.^(1/2);y2=x.^(1/3);y3=x.^(1/4);
>> plot(x,y1)
>> hold on
>> plot(x,y2)
>> hold on
>> plot(x,y3)


圖形一

(2) line指令和plot指令皆可以單一plot指令同時畫出三個圖形,但是兩者指令寫法略有不同。以line指令讓兩函數同時出現在一個指令時,要安排成行向量,每一行代表一條曲線:

line指令
x=1:1000;
y=line([x',x',x'],[x'.^(1/2),x'.^(1/3),x'.^(1/4)])

plot指令
x=1:1000;
plot(x,x.^(1/2),x,x.^(1/3),x,x.^(1/4))

而兩者畫出來的圖形會以三個不同顏色表示:藍→綠→紅(依先後順序排列)

圖形二

(3) line指令沒有plot指令之簡易設定功能,所以有些基本參數需依其參數名稱設定。如果要更改以line指令繪出多個圖形的顏色、線條寬度,或是線條形式,要用set呼叫的方式一一變更;而plot指令則可直接以單一指令指定曲線的形式。

line指令
>> x=1:10;
>> y1=x.^(1/2);y2=x.^(1/3);y3=x.^(1/4);
>> h1=line(x,y1);
>> h2=line(x,y2);
>> h3=line(x,y3);
>> set(h1,'color','r','linestyle','-','marker','o')
>> set(h2,'color','g','linestyle','-.','marker','x')
>> set(h3,'color','k','linestyle',':','marker','*')

plot指令
x=1:10;
plot(x,x.^(1/2),'r-o',x,x.^(1/3),'g-.x',x,x.^(1/4),'k:*')

圖形三

星期一, 11月 27, 2006

MATLAB課堂作業: 作業八

1. Find the multiplication result of two polynorminals, in which p=133x^5+122x^3+1, q=2x^4+100x^2+1.
ANS:
>> p=[133 0 122 0 0 1];
>> q=[2 0 100 0 1];
>> M=conv(p,q)

M =

266 0 13544 0 12333 2 122 100 0 1
故兩多項式乘積M= 266x9+13544x7+12333x5+2x4+122x3+100x2+1

2. A polynorminal is defined as f=100x^3+23x^2+x+45. Find the value f(x) if x is a magic matrix in the order of 5.
ANS:
>> x=magic(5);
>> f=[100 23 1 45];
>> v=polyval(f,x)

v =

498009 1395717 169 52725 342735
1228935 13125 35479 278967 415549
6817 22479 223645 809265 1075999
102355 176169 694267 936309 2955
135939 590715 1576945 939 74817

3. Using p and q defined in the item 1, find the quotient(商數) and residue(餘數) of p/q.
ANS:
>> p=[133 0 122 0 0 1];
>> q=[2 0 100 0 1];
>> [s,r]=deconv(p,q) % s代表商數,r代表餘數
s =

66.5000 0


r =

1.0e+003 *

0 0 -6.5280 0 -0.0665 0.0010

商數為66.5x,餘數為-6528x3-66.5x+1

4. Find the roots of p=0 and q=0, in which both p & q are defined in item 1.
ANS:
>> p=[133 0 122 0 0 1];
>> q=[2 0 100 0 1];
>> roots(p)
ans =

-0.0045 + 0.9578i
-0.0045 - 0.9578i
0.1039 + 0.1744i
0.1039 - 0.1744i
-0.1988

>> roots(q)
ans =

0 + 7.0704i
0 - 7.0704i
0 + 0.1000i
0 - 0.1000i

p共有5個根,q共有4個根

5. Fit a polynorminal curve to following data to an order of 3. x=[1:10]; y=[1210, 1866, 2301, 2564, 2724, 2881, 2879, 2915, 3010].
ANS:
>> x=[1:10];y=[1210, 1866, 2301, 2564, 2724, 2881, 2879, 2915, 3010];
>> polyfit(x,y,3)
??? Error using ==> polyfit
X and Y vectors must be the same size.

原本題目中y陣列中的元素值只有9個,與x陣列元素個數不符,因此執行指令時polyfit會有錯誤指令出現
若再加入一個自訂數值3178,則可算出最符合的多項式:
>> x=[1:10];y=[1210,1866,2301,2564,2724,2881,2879,2915,3010,3178];
>> polyfit(x,y,3)

ans =

1.0e+003 *

0.0065 -0.1372 1.0049 0.3410

多項式為6.5x3-137.2x2+100.49x+341

星期四, 11月 23, 2006

matlab課堂筆記 061124_第九章 繪圖

前言:
繪圖之一般程序可歸納如下:
(1)準備資料
(2)選擇適當圖形視窗及繪製區間
(3)呼叫基本繪圖函數
(4)選擇曲線型式及標記特性
(5)設定軸範圍及網線
(6)加註座標名稱、圖標及文字
(7)輸出圖形

繪圖指令 line & plot
繪線必須有起始點及終止點,但若有連續的點,採用矩陣表示時,則會自動連線,繪出整個圖形。
h(handle)=line(X,Y)
同時繪好幾條線時,可以以多行對應表示,每一行代表一線。
(如果用plot指令,需要同時出現兩個圖的話則繪圖指令中間需加上hold on指令)

* 如果要更改已經繪出圖形的顏色,則要用呼叫的方式 (y1為已繪出圖形的代碼handle)
set(y1,'Color','red')
改變線條寬度:
set(y2,'linewidth',6)
改變線條型式:
set(y3,'linestyle',':') (':'代表虛線)

* 以plot指令同時畫出三個圖形:
plot(x,x.^2,x,x.^3,x,x.^4) (畫出來的圖形就會以三個不同顏色表示:藍→綠→紅)
亦可直接在plot指令中指定線的形式:
plot(x,x.^2,'r-',x,x.^3,'go',x,x.^4,'k:')

* 若要在圖形中加上格線,可用grid on指令

* 線段指令亦可應用於三維的點線:
h=line(X,Y,Z)
h=line(X,Y,Z,'PropertyName',PropertyValue,...)


* 用line指令若要同時繪出線段的標記和連線,要畫兩次 (用plot指令則無此問題)

參數值的擷取與設定 get & set
get()set()為改變參數值可用設定指令,其格式如下:
set(sys,'Property1',Value1,'Property2',Value2,...)
or
a= get(h,'DefaultObjectTypePropertyName')


* 指令axis則有多重之應用,最主要是設定兩座標軸之範圍,其值是一個具有四個元素之列矩陣,前兩者分別為X軸之最小及最大值;後者為Y軸之最小與最大值。
axis([xmin xmax ymin ymax])
此外,axis之指令可以在後面加下面之中任一個參數:
axis auto/manual/tight/fill/ij/xy/equal/image/square/off/on

函數繪圖指令 ezplot & fplot
EZPLOT→ 以一般簡單函數進行繪圖。
ezplot(f) 是繪製函數 f(x) 之關係,其自變數x則需給予適當範圍。若在指令的參數項中未陳述範圍時[a,b],其預設值為 -2*pi < x < 2*pi。若要自行設定範圍,則可在f之後增加一個列矩陣參數,以表示其區間:
ezplot(f, [a,b])
p.s. ezplot3則為三維繪圖指令
ezplot3(x(t),y(t),z(t),[tmin tmax])
FPLOT→ 特定函數若具有M檔案之名稱,以 y = f(x)之型式,可以使用fplot 之指令進行繪圖。
在指定函數名稱時,可以使用字串指令eval,如:'sin(x)'、'[sin(x), cos(x)]',亦可使用函數手把呼叫。
fplot(function,limits,LineSpec)

併圖指令 Subplot
同一頁面要繪製許多圖時,可以使用subplot(m,n,p)這個併圖指令。
先宣告要畫幾個圖:(m,n)代表依m x n之矩陣型式分為m x n個小圖,p代表要畫的圖的位置。
要顯示3個圖的話,一樣要先宣告2x2的四分圖表示法,然後在subplot指令之最後項採用向量型式,如subplot(2,2,[1 2]) 表示上兩圖合併,而subplot(2,2, [1 3])表示左兩圖合併。
subplot只適用於2維圖形。

星期五, 11月 17, 2006

MATLAB課堂作業: 作業七

1. Write a program to input the name, age, sex and email address of a student and store the data as a structural format.

ANS:
程式指令:

% Q7_1.m
% input the name, age, sex and email address of a student
a=input('請輸入您的姓名:','s');
if isempty(a)
a='missing';
end
b=input('請輸入您的年齡:');
if isempty(b)
b='missing';
end
c=input('請輸入您的性別:(M/F)','s');
if isempty(c)
c='missing';
end
d=input('請輸入您的郵件地址:(ex. r94428004@ntu.edu.tw)','s');
if isempty(d)
d='missing';
end
student.Name=a;
student.Age=b;
student.Sex=c;
student.Email=d


執行結果:
(1)輸入資料齊全情況下
>> Q7_1
請輸入您的姓名:Susan
請輸入您的年齡:25
請輸入您的性別:(M/F)F
請輸入您的郵件地址:(ex. r94428004@ntu.edu.tw)r94428004@ntu.edu.tw

student =

Name: 'Susan'
Age: 25
Sex: 'F'
Email: 'r94428004@ntu.edu.tw'

(2)沒有輸入資料情況下則會輸出"missing"訊息
>> Q7_1
請輸入您的姓名:Susan
請輸入您的年齡:
請輸入您的性別:(M/F)F
請輸入您的郵件地址:(ex. r94428004@ntu.edu.tw)

student =

Name: 'Susan'
Age: 'missing'
Sex: 'F'
Email: 'missing'

2. Use the "menu" and "switch" commands incorporated in a program that enables a selection of Apple, Microsoft, IBM, Acer and Asus computers.

ANS:
程式指令:
% Q7_2.m
% a selection of Apple, Microsoft, IBM, Acer and Asus computers
com=menu('請選擇您想購買的電腦廠牌:','Apple','Microsoft','IBM','Acer','Asus');
switch com
case 1
fprintf('您選擇的電腦廠牌為 Apple\n')
case 2
fprintf('您選擇的電腦廠牌為 Microsoft\n')
case 3
fprintf('您選擇的電腦廠牌為 IBM\n')
case 4
fprintf('您選擇的電腦廠牌為 Acer\n')
case 5
fprintf('您選擇的電腦廠牌為 Asus\n')
end

執行結果:
在指令窗輸入Q7_2之後會出現一選擇視窗
選擇視窗

點選Acer之後則會在指令窗出現以下訊息:
您選擇的電腦廠牌為 Acer

3. Using fprintf and fscanf commands to store a magic matrix of order 5 in a file name magic.dat, and retrieve back to a set variable.

ANS:
magic(5)的矩陣如下
>> magic(5)

ans =

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

利用fprintf指令將magic(5)內容儲存於另一個檔案magicr94428004.dat
>> m=magic(5);
>> fid=fopen('magicr94428004.dat','w');
>> fprintf(fid,'%2.0f %2.0f %2.0f %2.0f %2.0f\n',m(1,:),m(2,:),m(3,:),m(4,:),m(5,:));
>> fclose(fid);
>> type magicr94428004.dat

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

輸出檔案結果如下:
輸出檔案

利用fscanf指令將檔案內容擷取於指令窗中
>> B=fopen('magicr94428004.dat');
>> A=fscanf(B,'%d',[5 5]);
>> fclose(fid);
>> A

A =

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

>> A'

ans =

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



4. Write a program to retrieve data from an excell sheet.

ANS:
首先先建立一個excell檔案,再將檔案內內容擷取出來

程式指令:
% Q7_4.m
F={'水果名稱' '產季' '產地' '價錢(斤)';'柳丁' '冬' '豐原' 20;'芒果' '夏' '玉井' 30;...
'草莓' '冬' '大湖' 50;'荔枝' '夏' '台南' 45;'蘋果' '一年四季' '梨山' 25};
xlswrite('D:\Physical Therapy\Graduate School\Class\碩二上\MATLAB Class\fruit.xls',F,'Sheet1');
[N,T,rawdata]=xlsread('D:\Physical Therapy\Graduate School\Class\碩二上\MATLAB Class\fruit.xls')

執行結果:

建立excel檔案與存放路徑
擷取結果
其中N為數值陣列,T為字串陣列,rawdata則包含所有資料內容(為一細胞陣列)

星期四, 11月 16, 2006

MATLAB課堂作業: 作業六

1. Run the function "freebody " in section 6.2 finding the height of the ball in terms of elapsed time in an interval of 0 to 30 seconds with an increment of 0.2 second.

ANS:
程式指令:

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(0:0.2:30)

ans =

1.0e+005 *

Columns 1 through 11

0 0.0002 0.0008 0.0018 0.0031 0.0049 0.0071 0.0096 0.0125 0.0159 0.0196

Columns 12 through 22

0.0237 0.0282 0.0331 0.0384 0.0441 0.0502 0.0566 0.0635 0.0708 0.0784 0.0864

Columns 23 through 33

0.0949 0.1037 0.1129 0.1225 0.1325 0.1429 0.1537 0.1648 0.1764 0.1884 0.2007

Columns 34 through 44

0.2134 0.2266 0.2401 0.2540 0.2683 0.2830 0.2981 0.3136 0.3295 0.3457 0.3624

Columns 45 through 55

0.3795 0.3969 0.4147 0.4330 0.4516 0.4706 0.4900 0.5098 0.5300 0.5506 0.5715

Columns 56 through 66

0.5929 0.6147 0.6368 0.6593 0.6823 0.7056 0.7293 0.7534 0.7779 0.8028 0.8281

Columns 67 through 77

0.8538 0.8798 0.9063 0.9332 0.9604 0.9880 1.0161 1.0445 1.0733 1.1025 1.1321

Columns 78 through 88

1.1621 1.1925 1.2232 1.2544 1.2860 1.3179 1.3502 1.3830 1.4161 1.4496 1.4835

Columns 89 through 99

1.5178 1.5525 1.5876 1.6231 1.6589 1.6952 1.7319 1.7689 1.8063 1.8442 1.8824

Columns 100 through 110

1.9210 1.9600 1.9994 2.0392 2.0794 2.1199 2.1609 2.2023 2.2440 2.2861 2.3287

Columns 111 through 121

2.3716 2.4149 2.4586 2.5027 2.5472 2.5921 2.6374 2.6830 2.7291 2.7756 2.8224

Columns 122 through 132

2.8696 2.9173 2.9653 3.0137 3.0625 3.1117 3.1613 3.2113 3.2616 3.3124 3.3636

Columns 133 through 143

3.4151 3.4670 3.5194 3.5721 3.6252 3.6787 3.7326 3.7869 3.8416 3.8967 3.9521

Columns 144 through 151

4.0080 4.0643 4.1209 4.1779 4.2354 4.2932 4.3514 4.4100

2. Draw two circles with radii of 3 and 5 cm each and their centers located at (0,0) and (10, 0) cm.

ANS:
因為欲求的圓心座標與半徑長已知,所以我以一個簡單的敘述指令執行

程式指令:
% drawcircles.m
% Draw circles with radii of r.
x1=0;
y1=0;
r1=3; % 第一個圓心座標為(0,0),半徑長為3的圓
theta=linspace(0,2*pi,60);
xx1=x1+r1*cos(theta);
yy1=y1+r1*sin(theta);
axis equal;
plot(xx1,yy1)
hold on; % 讓兩個圓同時出現在一張圖裡
x2=10;
y2=0;
r2=5; % 第二個圓心座標為(10,0),半徑長為5的圓
theta=linspace(0,2*pi,60);
xx2=x2+r2*cos(theta);
yy2=y2+r2*sin(theta);
axis equal;
plot(xx2,yy2)


執行結果:



3. Using an implicit function(匿名函數)to find the result for the formula y=x²/a²+y²/b², in which a=3, b=5, and y=[12 15 20] for x=[3 4 5]. Draw the graphs.

ANS:
>> a=3;b=5;g=@(x,y) x.^2/a^2+y.^2/b^2;
>> g([3 4 5],[12 15 20])

ans =

6.7600 10.7778 18.7778

4. A plate ring has a thickness of 2 cm, with inner and outer diameters of 3 cm and 5 cm, respectively. Using the above parameters as inputs of a function to find the area of its upper face and side wall as well as the volume.

ANS:
程式指令:
function [area volume]=ring(Do,Di,height)
% Find the area and volume of a plate ring
% Do,Di:outside & inside diameters
area=abs(Do.^2-Di.^2)*pi;
volume=abs(Do.^2-Di.^2).*height*pi;

執行結果:
>> [area volume]=ring(5,3,2)

area =

50.2655


volume =

100.5310

5. Rewrite the function "nest_fun" in the text by moving the nested function polyx(xx) outside of the main function so that the function will still produce the same results.

ANS:
原有之程式指令為:
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

執行結果:
>> nest_fun(2:10,[1 2 4;2 4 8])

ans =

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

若要將polyx(xx)移至主函數程式之外,則需要對A、B、C變數加上全域性變數(global)之宣告:
function [rr_array]=nest_fun2(x,a)
%function to find sets of polynormials.
% a: set of constants, [A B C]
% x: variables in array
% Example: rr=nest_fun2(2:10,[1 2 4;2 4 8])
n=size(a);
for i=1:n
global A B C
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
end
function r=polyx(x)
global A B C
r=A.*x.^2 + B.*x +C;
end

執行結果:
>> nest_fun2(2:10,[1 2 4;2 4 8])

ans =

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

與原程式之結果相同

6. rewrite the function "pillar" to make sure it will run smoothly even the inputs are not sufficient for the function to run.

ANS:
原有程式只適用於完全沒有輸入數值時仍可有結果輸出,
但若輸入數值為一個或兩個時,則執行原程式時會有error發生,
故須將程式改為:
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;
elseif nargin<2,
height=1;
Di=0;
elseif nargin<3,
height=1;
end
volume=abs(Do.^2-Di.^2).*height*pi/4;

執行結果:
>> pillar

ans =

0.7854

>> pillar(2)

ans =

3.1416

>> pillar(1,2)

ans =

2.3562

>> pillar(1,2,2)

ans =

4.7124

如此,不論輸入數值是否小於3個都可以順利執行程式

星期六, 11月 11, 2006

MATLAB課堂作業: 作業五

1. Create a 3X3X2 matrix A which contains a matrix of magic(3) and another rand(3,3)*10.

ANS:
>> A(:,:,1)=magic(3); % 第一頁為magic(3)矩陣
>> A(:,:,2)=rand(3,3)*10 % 第二頁為rand(3,3)*10矩陣

A(:,:,1) =

8 1 6
3 5 7
4 9 2


A(:,:,2) =

4.4470 9.2181 4.0571
6.1543 7.3821 9.3547
7.9194 1.7627 9.1690

2. Assume that matrices A=[45 89 99; 12 34 55], B=[15 25 45; 65 50 30]. Find matrices that join A & B both in vertical and horizontal directions. Also, find a two-page matrix with A & B stored in separate pages. Use commands horizcat and verticat to check the results with those from above.

ANS:
利用cat指令將兩矩陣合併
>> A=[45 89 99; 12 34 55];
>> B=[15 25 45; 65 50 30];
>> Cvertical=cat(1,A,B) % dim=1時,表示將A、B矩陣垂直疊加

Cvertical =

45 89 99
12 34 55
15 25 45
65 50 30

>> C1=[A;B]

C1 =

45 89 99
12 34 55
15 25 45
65 50 30

兩者結果相同

>> Chorrizontal=cat(2,A,B) % dim=2時,表示將A、B矩陣水平疊加

Chorrizontal =

45 89 99 15 25 45
12 34 55 65 50 30

>> C2=[A,B]

C2 =

45 89 99 15 25 45
12 34 55 65 50 30

兩者結果相同

>> Cpage=cat(3,A,B) % dim=3時,表示將A、B矩陣以分置於不同頁方式疊加

Cpage(:,:,1) =

45 89 99
12 34 55


Cpage(:,:,2) =

15 25 45
65 50 30

3. Find a matrix M which contains A as the first and third pages and B the second and fourth pages.

ANS:
>> A=[45 89 99; 12 34 55];
>> B=[15 25 45; 65 50 30];
>> M=cat(3,A,B,A,B)

M(:,:,1) =

45 89 99
12 34 55


M(:,:,2) =

15 25 45
65 50 30


M(:,:,3) =

45 89 99
12 34 55


M(:,:,4) =

15 25 45
65 50 30

4. Construct a 2X2 cell array that contains 'Eric' [90 100]; 'Peter' [35 60]; 'Jan' [77 67].

ANS:
>> D=cell(3,2);
>> D={'Eric' [90 100]; 'Peter' [35 60]; 'Jan' [77 67]}

D =

'Eric' [1x2 double]
'Peter' [1x2 double]
'Jan' [1x2 double]

5. Construct a structural Array that contains the following data:
Name: 'Philip', 'Peter','Roberts', 'Roe'
Age: 35, 45, 55, 60
Salary: 50,000 40,000 80,000 120,000


ANS:
>> Database(1).Name='Philip';Database(1).Age=35;Database(1).Salary=50000;
>> Database(2).Name='Peter';Database(2).Age=45;Database(2).Salary=40000;
>> Database(3).Name='Roberts';Database(3).Age=55;Database(3).Salary=80000;
>> Database(4).Name='Roe';Database(4).Age=60;Database(4).Salary=120000;
>> Database % 顯示整個結構性矩陣

Database =

1x4 struct array with fields:
Name
Age
Salary

亦可顯示結構性矩陣內各筆資料內容
>> Database(1) % 第一筆

ans =

Name: 'Philip'
Age: 35
Salary: 50000

>> Database(2) % 第二筆

ans =

Name: 'Peter'
Age: 45
Salary: 40000

>> Database(3) % 第三筆

ans =

Name: 'Roberts'
Age: 55
Salary: 80000

>> Database(4) % 第四筆

ans =

Name: 'Roe'
Age: 60
Salary: 120000

另外,以struct指令亦可得到相同的結果
>>Database2=struct('Name',{'Philip','Peter','Roberts','Roe'},'Age',{35,45,55,60},'Salary',{50000,40000,80000,120000})

Database2 =

1x4 struct array with fields:
Name
Age
Salary

>> Database2(1)

ans =

Name: 'Philip'
Age: 35
Salary: 50000

>> Database2(2)

ans =

Name: 'Peter'
Age: 45
Salary: 40000

>> Database2(3)

ans =

Name: 'Roberts'
Age: 55
Salary: 80000

>> Database2(4)

ans =

Name: 'Roe'
Age: 60
Salary: 120000

星期五, 11月 10, 2006

MATLAB課堂作業: 作業3.3

五月份的台北天氣已經開始轉熱,其整月的日平均溫度如下:
T=[25 26 24 23, 18 17 20 30 29, 25 19 19 21 20, 23 24 23 25 23,…
27 28 24 26 27, 23 25 24 21 19, 20]C
1. 試撰寫一程式,以上述之溫度資料為輸入值,求其月平均、最高及最低溫度。
2. 試撰寫另一程式,以交談式輸入之方式,按月份輸入每日之溫度,(必須分清楚每月應有之輸入日數)然後呼叫第一項之程式進行分析。


ANS:
1.
程式指令
function [Mean Max Min]=temp(T)
% 計算月平均、最高及最低溫度
Mean=mean(T); % 平均溫度
Max=max(T); % 最高溫度
Min=min(T); % 最低溫度

執行結果
>> T=[25 26 24 23 18 17 20 30 29 25 19 19 21 20 23 24 23 25 23 27 28 24 26 27 23 25 24 21 19 20];
>> [Mean Max Min]=temp(T)

Mean =

23.2667


Max =

30


Min =

17


2.
程式指令

% month_temp.m
% 以交談方式,按月份輸入每日之溫度,呼叫另一程式進行溫度分析
month=input('請輸入月份 (1-12):');
T=input('請輸入每日溫度 ex.[24 28 30]:');
switch month
case {1,3,5,7,8,10,12}
if numel(T)==31
[Mean Max Min]=temp(T)
else
input('輸入溫度之日數不符,此月份應輸入31天之溫度!\n')
end
case 2
if numel(T)==28
[Mean Max Min]=temp(T)
else
input('輸入溫度之日數不符,此月份應輸入28天之溫度!\n')
end
case {4,6,9,11}
if numel(T)==30
[Mean Max Min]=temp(T)
else
input('輸入溫度之日數不符,此月份應輸入30天之溫度!\n')
end
end


執行結果
& 情境一 (當輸入溫度之日數與當月天數相符時,則可算出當月平均、最高及最低溫度):
>> month_temp
請輸入月份 (1-12):4
請輸入每日溫度 ex.[24 28 30]:[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30]

Mean =

15.5000


Max =

30


Min =

1

& 情境二 (當輸入溫度之日數與當月天數不符時,則會顯示警示):
>> month_temp
請輸入月份 (1-12):12
請輸入每日溫度 ex.[24 28 30]:[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30]
輸入溫度之日數不符,此月份應輸入31天之溫度!

MATLAB課堂作業: 作業4.3

氣體在特定容器下,其壓力、體積與溫度間之關係如下:P=nRT/V。這是波義耳的氣體公式,但萬得爾(van der Vaals)公式則加入分子間之引力因素。其關係如下:P=nRT/(V-nb)-an²/V²。其中,a,b分別為氣體之特定常數。理想氣體在n=1摩爾數下,溫度在0C(273.2K)下,氣壓為一大氣壓時,具有體積為22.4公升,其R值為0.0826。對於氯氣而言,上述之a,b值分別為a=6.49,b=0.0562。試就上兩公式計算其壓力差。

ANS:
程式指令

% Q4_3.m
% 分別以波以耳氣體公式與萬得爾公式計算一莫耳的理想氣體和氯氣的壓力差
% P=nRT/V 波以耳氣體公式 P=nRT/(V-nb)-an2/V2 萬得爾公式
% n=1,R=0.0826,T=273.2,V=22.4,a=6.49,b=0.0562
P1=1*0.0826*273.2/22.4;
P2=1*0.0826*273.2/(22.4-1*0.0562)-6.49*1^2/22.4^2;
P=P1-P2;
absp=abs(P);
if P>0
fprintf('波以耳公式算出之理想氣體壓力值較萬得爾公式算出之氯氣壓力值大 %0.4f\n',absp)
elseif P<0
fprintf('波以耳公式算出之理想氣體壓力值較萬得爾公式算出之氯氣壓力值小 %0.4f\n',absp)
else
fprintf('波以耳公式算出之理想氣體壓力值與萬得爾公式算出之氯氣壓力壓力值相等\n')
end


執行結果
>> Q4_3
波以耳公式算出之理想氣體壓力值較萬得爾公式算出之氯氣壓力值大 0.0104

星期四, 11月 09, 2006

HTML語法:表格製作

表格製作




姓名學號成績
河馬1234567200


註解:
tr 行
th 欄
align="center" 置中

星期五, 11月 03, 2006

主函數與次函數的應用:畫圓

畫圓的函數檔
以主函數與次函數構成

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


註1:執行時在command window中輸入
>>drawcircles(m:n)
即可畫出圖形
再輸入
>>hold off
>>figure(2)
>>drawcircles(m:n)
則會在另一個視窗中(figure2)畫出新圖形
因為程式中以隨機出現圓心位置和半徑,因此每次畫出來的圖都不一樣
註2: 在指令中加入rand('state',0)宣告,則每次畫出來的圖都會一樣

星期四, 11月 02, 2006

MATLAB課堂作業: 作業四

1. Two vectors A=3i +5i +10k and B=5i -6j +2k are both passing through an origin. Find the unit vector that is perpendicular to both vectors A & B.

ANS:
可以用外積的公式算出與兩向量垂直之向量
若A=a1i+a2j+a3k且B=b1i+b2j+b3k
則垂直向量C=A×B=(a2*b3-a3*b2)i+(a3*b1-a1*b3)j+(a1*b2-a2*b1)k
可以cross指令計算得出:
>> A=[3 5 10];B=[5 -6 2];
>> C=cross(A,B)

C =

70 44 -43

故向量C=70i+44j-43k

2. A=randn(10,10). Find mean(A), median(A,2), std(A), std(A,1,2). Explain the results.

ANS:
>> A=randn(10,10)

A =

-1.1878 -1.1859 0.1286 0.8057 -0.3306 -0.1199 0.4694 1.0184 -0.4650 1.5532
-2.2023 -1.0559 0.6565 0.2316 -0.8436 -0.0653 -0.9036 -1.5804 0.3710 0.7079
0.9863 1.4725 -1.1678 -0.9898 0.4978 0.4853 0.0359 -0.0787 0.7283 1.9574
-0.5186 0.0557 -0.4606 1.3396 1.4885 -0.5955 -0.6275 -0.6817 2.1122 0.5045
0.3274 -1.2173 -0.2624 0.2895 -0.5465 -0.1497 0.5354 -1.0246 -1.3573 1.8645
0.2341 -0.0412 -1.2132 1.4789 -0.8468 -0.4348 0.5529 -1.2344 -1.0226 -0.3398
0.0215 -1.1283 -1.3194 1.1380 -0.2463 -0.0793 -0.2037 0.2888 1.0378 -1.1398
-1.0039 -1.3493 0.9312 -0.6841 0.6630 1.5352 -2.0543 -0.4293 -0.3898 -0.2111
-0.9471 -0.2611 0.0112 -1.2919 -0.8542 -0.6065 0.1326 0.0558 -1.3813 1.1902
-0.3744 0.9535 -0.6451 -0.0729 -1.2013 -1.3474 1.5929 -0.3679 0.3155 -1.1162

>> mean(A) %結果顯示每一行的平均值

ans =

-0.4665 -0.3757 -0.3341 0.2245 -0.2220 -0.1378 -0.0470 -0.4034 -0.0051 0.4971

>> median(A,2) %結果顯示每一列的中間值

ans =

0.0044
-0.4545
0.4915
-0.2024
-0.2061
-0.3873
-0.1415
-0.4096
-0.4338
-0.3712

>> std(A) %結果顯示每一行的標準差,其計算式為以(n-1)為底的公式做計算

ans =

0.9150 0.9908 0.7798 0.9822 0.8485 0.7584 0.9899 0.7716 1.1289 1.1608

>> std(A,1,2) %結果顯示每一列的標準差,其計算式為以n為底的公式做計算

ans =

0.8629
0.9447
0.9378
0.9868
0.9209
0.8240
0.8131
1.0327
0.7426
0.9075

3. Prove that exp(i*theta)=cos(theta)+i*sin(x) using matlab commands.

ANS:
>> theta=[30 40]*180/pi; % 以角度為30度和40度做驗證,先換算為弧度
>> a=exp(i*theta);
>> b=cos(theta)+i*sin(theta);
>> isequal(a,b)

ans =

1

故exp(i*theta)=cos(theta)+i*sin(theta)得證
p.s. 兩個角度算出來的值分別為(-0.9122 - 0.4098i)和(0.0393 - 0.9992i)

4. Let R=eye(3)*5+4*ones(3)i. Find abs(R), angle(R), real(R) and imag(R).

ANS:
>> R=eye(3).*5+4*ones(3)*i

R =

5.0000 + 4.0000i 0 + 4.0000i 0 + 4.0000i
0 + 4.0000i 5.0000 + 4.0000i 0 + 4.0000i
0 + 4.0000i 0 + 4.0000i 5.0000 + 4.0000i

>> abs(R)

ans =

6.4031 4.0000 4.0000
4.0000 6.4031 4.0000
4.0000 4.0000 6.4031

>> angle(R) % 弧度

ans =

0.6747 1.5708 1.5708
1.5708 0.6747 1.5708
1.5708 1.5708 0.6747

>> angle(R)*180/pi % 角度

ans =

38.6598 90.0000 90.0000
90.0000 38.6598 90.0000
90.0000 90.0000 38.6598

>> real(R) % the real part of the elements

ans =

5 0 0
0 5 0
0 0 5

>> imag(R) % the imaginary part of the elements

ans =

4 4 4
4 4 4
4 4 4

MATLAB課堂作業: 作業三

1. for i=1:50,
A(i)=i^2;
end;
A=reshape(A,5,10)
Find A. Is there any better way to find the same result?


ANS:
指令亦可改成
>> A=[1:50].^2;
>> A=reshape(A,5,10);
>> find(A)

ans =

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
如果只是要得到相同的結果,則亦可省略reshape的指令
>> A=[1:50].^2;
>> find(A)
但是單就如此,其結果顯示的矩陣排列方式和原本的結果不同
因此指令須改為
>> A=[1:50].^2;
>> find(A)’

2. Let A=magic(5). Replace the elements that are larger than 10 with 'NaN'

ANS:
>>clear all
>> 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

>> A(A>10)=NaN

A =

NaN NaN 1 8 NaN
NaN 5 7 NaN NaN
4 6 NaN NaN NaN
10 NaN NaN NaN 3
NaN NaN NaN 2 9

3. Replace the prime elements in Prob 2 with 0.

ANS:
本來以為題意是要將已經做過取代動作的矩陣中的質數以0取代,但是因為新矩陣中有非整數出現,因此無法執行isprime的指令找出質數。若是先將原本A矩陣中的質數代為0,再將新矩陣中>10的數值以NaN取代,其求出的解亦與原意不同,故不適用。若只是單純把原矩陣A中的質數代為0,其指令和結果如下:
>> clear
>> 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

>> A(isprime(A))=0

A =

0 24 1 8 15
0 0 0 14 16
4 6 0 20 22
10 12 0 21 0
0 18 25 0 9

4. The function draw_number is shown in sec 3.3. Please make a comparison of results when no_of_draw=[100:100:1000]?

ANS:
將no_of_draw=[100:100:1000]代入sec 3.3程式後執行所得之結果為
>> draw_number([100:100:1000])

ans =

17 18 20 21 24

由結果可知其抽籤之總次數為100次
和no_of_draw=100代入有相似之結果
>> draw_number(100)

ans =

17 24 20 20 19

其原因為:
在程式中的指令是在while n<=no_of_draw的前提下重複執行
而no_of_draw=[100:100:1000]相當於[100 200 300 400 500 600 700 800 900 1000]之陣列
當n<=100時,陣列中所有的條件皆成立,
>> 100<=[100:100:1000]

ans =

1 1 1 1 1 1 1 1 1 1
但是當n>100時,則陣列中至少有其中一個元素不成立
>> 101<=[100:100:1000]

ans =

0 1 1 1 1 1 1 1 1 1
因此以no_of_draw=[100:100:1000]代入程式中計算時,只會運算到第100次

5. Please use while...end statement to find out the chances to obtain a preset two-digit number which will be similar to the number set drawn from function rand().

ANS:
程式指令:
% chances_draw_number.m
% use while...end statement to find out the chances to obtain a preset
% two-digit number which will be similar to the number set drawn from
% function rand()

num=input('請輸入一個十位數字:');
no_of_draw=input('請輸入欲計算輸入數字與隨機號碼相同之次數:');
n=0;
for i=1:no_of_draw
while 1
ball=fix(rand*100);
if ball~=num
n=n+1;
else
break;end
chance(i)=1/n;
end
allchance=sum(chance(i))/no_of_draw*100;
end
fprintf('輸入數字 %d,與隨機號碼相同次數達 %d次的機率為 %0.3f%%\n',num,no_of_draw,allchance)

執行結果:
>> chances_draw_number
請輸入一個十位數字:12
請輸入欲計算輸入數字與隨機號碼相同之次數:5
輸入數字 12,與隨機號碼相同次數達 5次的機率為 0.059%

林書如:您的體重過重嗎?

前言:
BMI指數可以反映出體重是否處於理想狀態。過去已有許多文獻指出BMI指數過高往往是造成糖尿病、高血壓,甚至是心血管疾病的危險因子。因此,適度地監控自己的BMI指數是否處於正常值內是很重要的。

問題:
試寫出一個函數,當輸入體重和身高的數值後,能計算出其身體質量指數BMI值。並且依據其BMI值判斷其是否有體重過重或過輕的問題。(18.5<=BMI<24時顯示"體重標準",BMI>24時顯示"體重過重",BMI<18.5時顯示"體重過輕")

指令:

一開始我寫了一個比較簡單的函數檔(function file)

function [A]=BMI(w,h) 
% BMI(Body Mass Index) calculation
A=w/h^2; % A means BMI
if 18.5<=A & A<24
fprintf('BMI=%0.2f 體重標準\n',A)
elseif A<18.5
fprintf('BMI=%0.2f 體重過輕\n',A)
else
fprintf('BMI=%0.2f 體重過重\n',A)
end


模擬輸出結果:

>> BMI(60,1.5)
BMI=26.67 體重過重

ans =

26.6667

結果顯示,一位體重60公斤,身高150公分的人,其BMI指數為26.67,有體重過重的傾向

==============================================

但是,為了能夠加上交談的功能,讓使用者更能清楚地了解應該輸入的內容
例如,輸入之體重應以公斤為單位,輸入之身高應以公尺為單位
因此我又改寫成另一個敘述檔(Script file),如下

% BMI2.m
% BMI(Body Mass Index) calculation
a=input('Please keyin your body weight? (kg)');
if isnumeric(a)
b=input('Please keyin your body height? (m)');
B=a/b^2; % B means BMI
if 18.5<=B & B<24
fprintf('BMI=%0.2f 體重標準\n',B)
elseif B<18.5
fprintf('BMI=%0.2f 體重過輕\n',B)
else
fprintf('BMI=%0.2f 體重過重\n',B)
end
end


模擬輸出結果:

>> BMI2
Please keyin your body weight? (kg)60
Please keyin your body height? (m)1.5
BMI=26.67 體重過重

結果同上,一位體重60公斤,身高150公分的人,其BMI指數為26.67,有體重過重的傾向

討論:

我是以簡單的if-elseif指令來撰寫這兩個小程式
程式中如果還有可以補強的地方(例如發現程式的bias,或是其實還有更簡單的寫法),還麻煩大家提供意見 :p

星期四, 10月 26, 2006

html語法: 上標與下標

上標
a2+b2=c2
語法 : a2+b2=c2

下標
a2+b2=c2
語法 : a2+b2=c2

星期五, 10月 13, 2006

Auditory blog

網址 podcast.com

星期四, 10月 05, 2006

MATLAB課堂作業: 作業二

1. During SARS period, people who passed the airport gate should have their temperatures measured. Suppose that B is the recorded body temperatures in Fahrenheit and B=[98 100 102 104 98.4 98.2 98.5 101 102 99.5]. Please workout a program to do: a) Converse the unit of B into Centigrade b)Decide who got the temperature? Design a matrix D to show that "1" for the one who has a fever (>38C) and "0" for the normal.

ANS:
>> % Quesiotn 1
>>% 題旨:將華氏溫度轉換為攝氏溫度,並且定義出溫度大於38度c者為1
>> B=[98 100 102 104 98.4 98.2 98.5 101 102 99.5]; % 所有被測者的華氏體溫
>> c=(B-32)*5/9 % 將華氏溫度換算為攝氏溫度

c =

Columns 1 through 8

36.6667 37.7778 38.8889 40.0000 36.8889 36.7778 36.9444 38.3333

Columns 9 through 10

38.8889 37.5000

>> D=c>38 % 將攝氏溫度大於38度c者定義為1,得一新矩陣D

D =

0 0 1 1 0 0 0 1 1 0
>> % 執行結果: 換算後的攝氏溫度為[36.6667 37.7778 38.8889 40.0000 36.8889 36.7778 36.9444 38.3333 38.8889 37.5000]; 定義體溫大於38度c的matrix D=[0 0 1 1 0 0 0 1 1 0]
>> % 討論: 定義出高體溫者之後,若現實上受測者有編號,則可以利用尋找索引值的方式輕易找出高體溫者為誰

2. There are two matrices A and B, in which A=[ -10 8 6 4 -5 20] and B=[2 8 5 10 -6 3]. Find D using both Matlab commands .(a) D= (A>B) (b) D= (A>5) (f) D=A+B.

ANS:
>> % Quesiotn 2
>> % 題旨:比較矩陣A與矩陣B
>> A=[-10 8 6 4 -5 20];
>> B=[2 8 5 10 -6 3];
>> D=A>B % 矩陣A與矩陣B一對一比較,A>B者定義為1,其餘為0

D =

0 0 1 0 1 1

>> D=A>5 % 矩陣A中數值>5者定義為1,其餘為0

D =

0 1 1 0 0 1

>> D=A+B % 將矩陣A與矩陣B相加

D =

-8 16 11 14 -11 23

>> % 執行結果: (1) D=[0 0 1 0 1 1] (2) D= [0 1 1 0 0 1] (3) D=[-8 16 11 14 -11 23]
>> % 討論:矩陣A與矩陣B排列相同,因此兩者可以直接做比較與計算

3. Suppose that x=[10 20 30] and y=[1 4 6], Calculate the following combinations using Matlab commands:(a) A= 3x +y (b) B=5y/x (c) C=4x2y (d) D=sin(x)cos(y) (e)5x sin(2y)

ANS:
>> % Quesiotn 3
>> % 題旨:進行矩陣x與矩陣y之計算
>> x=[10 20 30];
>> y=[1 4 6];
>> A=3*x+y

A =

31 64 96

>> B=5*y./x

B =

0.5000 1.0000 1.0000

>> C=4*2*y

C =

8 32 48

>> D=sin(x).*cos(y)

D =

-0.2939 -0.5967 -0.9487

>> E=5*x.*sin(2*y)

E =

45.4649 98.9358 -80.4859
>> % 執行結果: A= [31 64 96], B= [0.5000 1.0000 1.0000], C=[8 32 48], D= [-0.2939 -0.5967 -0.9487], E= [45.4649 98.9358 -80.4859]
>> % 討論:重點在於元素對元素的計算符號是在矩陣的計算符號之前加上"."

4. A set of resistors are expressed in an array as R=[ 10 30 200 400]ohms. Write a M-file program to calculate the equivalent resistance if (a) all the resistors are in series (b) all of the resistors are in parallel.

ANS:
% Quesiotn 3
% 題旨:撰寫一個計算多個電阻分別串聯與並聯後之總電阻的M-file
R=[10 30 200 400];
R1=sum(R) % 串聯後之總電阻R1為四個電阻之總合 R1=r1+r2+r3+r4
R2=1/sum(1./R) % 並聯後之總電阻R2=1/(1/r1+1/r2+1/r3+1/r4)
>>% 在command window下執行M-file
>>% 執行結果: 串聯後之總電阻R1為640歐姆; 並聯後之總電阻R2為7.1006歐姆
>>% 討論:計算R矩陣中各個element的倒數不能單純以1/R指令執行,而需在/之前加上.才可順利算出數值

5. Suppose matrix A=1:24, please reform A into matrices B in size 3x8; C in size 6x4; D in size 2x12.

ANS:
>>% Question 5
>>% 將矩陣A(1x24)轉置為B (3x8), C(6x4), D(2x12)三種不同排列方式之矩陣
>> A=[1:24];
>> B=reshape(A,3,8)

B =

1 4 7 10 13 16 19 22
2 5 8 11 14 17 20 23
3 6 9 12 15 18 21 24

>> C=reshape(A,6,4)

C =

1 7 13 19
2 8 14 20
3 9 15 21
4 10 16 22
5 11 17 23
6 12 18 24

>> D=reshape(A,2,12)

D =

1 3 5 7 9 11 13 15 17 19 21 23
2 4 6 8 10 12 14 16 18 20 22 24
>>% 執行結果:如上所示之B、C、D三矩陣
>>% 討論:使用reshape(A,m,n)方式將陣列A的元素由上到下,由左到右的次序重新排列成一個新的m×n的矩陣。使用reshape函數時,重排後矩陣元素的個數必須與原矩陣相同。

星期六, 9月 30, 2006

MATLAB課堂作業: 作業一

1. The unit of POWER can be expressed as W, kW, hp, etc. Their equivalent expressions are J/s, kJ/s, in which J is joule to represent work and 1 J=0.101972kgm (You may find unit conversion table at http://www.digitaldutch.com/unitconverter/). Now, by using such relations find out the following condition when a farm tractor pulls a farm equipment working on the field at 20m/min. Suppose the draft force is meaured as 200kgf. Find out the power of the tractor.
(試利用上述的關係,計算下面的情況:設若一台曳引機(Tractor)拉引某一農具,其拉力測得200公斤,試計算此需要多少功率(動力、馬力或kW)才能該農具在一分鐘內拉動20m?)


ANS:
>> % Question 1
>> % To find out the power of the tractor
>> f=200; % 曳引機(tractor)拉力f為200公斤
>> d=20; % 向前拉動20m
>> J=0.101972; % 單位換算1 J=0.101972kgm
>> t=60; % 需時60秒
>> p=f*d/J/t % 功率p=f*d/J/t

p =

653.7742

>> % 執行結果:所需之功率為653.7742瓦特

2. A circuit with three resistors rated 1, 20,200 ohms in series is connected to a voltage source of 110V. Find the line current and power consumption.
(有一電路,以1、20及200歐姆之電阻器相串聯,試計算當此電路兩端加於110V電源時,其流經之電流及所耗之功率。)


ANS:
>> % Question 2
>> % Find the line current and power consumption
>> r=1+20+200; % 三電阻器串聯後之總電阻r
>> i=110/r; % 110伏特之電源除以總電阻得流經之電流i
>> p=i*110; % 電流與電壓之乘積得消耗功率
>> i

i =

0.4977

>> p

p =

54.7511
>> % 執行結果: 流經之電流為0.4977安培及所耗之功率為54.7511瓦特
>> % 討論: 第一、二題透過已有之物理基礎,了解消耗功率計算之方式,再加以MATLAB語法逐步計算求得解答,並非難事。

3. Calculate the area of a trangle with three sides in 650cm, 428cm, 282cm, repectively.
(一個三角形三邊分別為650、428、282cm,試計算其面積為多少平方公分。)


ANS:
>> % Question 3
>> % To calculate the area of a trangle
>> a=650;
>> b=428;
>> c=282; % 三角形三邊長分別為650、428、282公分
>> r=(a+b+c)/2;
>> area=(r*(r-a)*(r-b)*(r-c))^(1/2) % 以海龍公式計算三角形面積

area =

4.5233e+004

>> % 執行結果: 三角形面積為(4.5233e+004)平方公分

4. A hog house was designed the way that heat dissipations from lighting & equipment 25,000BTU/hr, heat flux through wall & cellings 52,000BTU/hr , heat flux through aerations 84,300BTU/hr. Suppose each hog will dissipate 600BTU/hr, how many hogs this house can raise?
(在一個規模的現代化豬舍中,其照明及設備操作所產生的熱量為25,000BTU;設經過牆壁及屋頂之流出之熱流量每小時為52,000BTU(英熱單位),利用風扇散熱的熱流量每小時為84,300BTU。然而飼養的豬隻中,每隻每小時產生600BTU之熱。試求此豬場要養多少隻豬才是合適的情況?)


ANS:
>> % Quesiotn 4
>> % How many hogs this house can raise?
>> a=52000+84300; % 每小時總散熱量為a
>> h=(a-25000)/600; % 每小時總散熱量a減去照明及設備操作所產生的熱量為所有豬隻能產生的最大熱量
>> % 所有豬隻能產生的最大熱量除以每隻豬每小時產熱後得一數值h
>> floor(h) % 取出小於或等於數值h的最大數目即為合適的豬隻量,以floor(x)指令執行

ans =

185
>> % 執行結果: 豬舍中養185隻豬最合適

5. The flow rate of a rectangular water gate can be calculated from the following formula: Q=3.33x[L-0.2H]^(3/2), in which L: width of gate in ft, H:water head in ft and Q: the flow rate in cu. ft./s. Find the flow rate Q if L=4 ft and the wate head is 0.9 ft.
(一個矩形堰之流量率通常用下列公式計計:Q=3.33x[L-0.2H]^(3/2),其中Q為流量(立方呎/秒);L為堰開口寬度(呎);H為堰口溢流水頭高度(呎)。若堰口寬度為4呎,堰口溢流水頭高度為0.9呎,問其流量每秒為多少立方呎?)


ANS:
>> % Question 5
>> % Find the flow rate Q
>> L=4; % 堰口寬度L為4呎
>> H=0.9; % 堰口溢流水頭高度H為0.9呎
>> Q=3.33*(L-0.2*H)^(3/2) % 流量為Q

Q =

24.8622
>> % 討論: 一剛開始看到已有公式,不加多想變複製至MATLAB視窗中執行,出現error訊息後才驚覺須以MATLAB語法重新闡述方可順利計算求解,程式語法之撰寫實在需要多加謹慎小心