以平行一級反應為例

Download Report

Transcript 以平行一級反應為例 

使用 MatLab 解微分方程
1
使用 MatLab 解微分方程式的基本流程
與相關基本知識
1. 首先撰寫描述該微分方程式的M檔案。 M檔案就像是一個小程式，
裡面可以是欲執行的命令、函數、運算式或程式語句。
2. 將該M檔案與初始條件值，載入ode求解函數中，透過求解器求出微
分方程的解。
3. 最後經過ode求解函數運算後，將傳回指定的時間區間內的解，以一
個行向量儲存。
4. 行向量是一個列數恆為1的矩陣，矩陣在MatLab中表示為(行,列) ，
例如: (1,1)、(2,1) …
(1,1)代表”1 x 1”的矩陣，可以儲存一個變數值；
(2,1)代表”2 x 1”的矩陣，可以儲存兩個變數值；以此類推。
5. 因此要解微分方程時，我們必須先撰寫好兩個M檔案，一個是
用來描述該微分方程式，一個則是用來進行初始條件值的設定。
2
Part I: M檔案的撰寫
(以一級反應為例)
一級反應的微分式為:
• Step 1:開啟MatLab 程式
點擊桌面上的MatLab圖示，執行MatLab。
(MatLab剛開啟時較慢，請耐心等候。)
d[A]/dt = -k*[A]
• Step 2:待紅色圈圈處的>>及Ready出現，即為開啟完成。
3
• Step 3:開啟 Editor 視窗
點選左上角 File → New → Function，會出現一個 Editor 視窗。
4
• Step 4:編輯描寫微分方程的M檔案
在Editor視窗中輸入以下指令:
function da=firstorder(t,a,k)
da=zeros(1,1);
da(1)=-k*a(1);
end
[這些指令下一頁會進行
詳細的說明]
5
function da=firstorder(t,a,k)
da=zeros(1,1);
da(1)=-k*a(1);
end
一級反應的微分式為:
d[A]/dt=-k*[A]
1. 檔案的內容中，第一行必須遵循以下格式:
function d輸出參數=函數名稱(輸入參數)
function
da
= firstorder
(t,a,k)
這裡 a 是輸出參數。
而其中函數名稱可自訂 ，但不可與變數名稱相同，且第一個
字元必須為字母，字母間不可空格。
這裡 輸入參數共有三個，分別為 t, a & k 。
6
function da=firstorder(t,a,k)
da=zeros(1,1);
da(1)=-k*a(1);
end
一級反應的微分式為:
d[A]/dt=-k*[A]
2. 第二行是將輸出參數“da “ 透過“zeros(1,1)” 將其設定為一個
1x1其值為零的矩陣。這個1x1的矩陣是要用來存放[A]值。所
以a(1)所代表的就是[A]。
3. 第三行則就用來描述一級反應的微分式， “ * ”號代表數學上的
乘號。
4. 第四行的end跟第一行開頭的function是互相對應的，代表該函
數運算的結束，但不能漏掉。
※在MatLab中，於程式行末尾加上分號，會使該行程式執行結果不顯示。
7
• Step 5:將M檔案存檔
點選左上角 File → Save As，儲存檔案，檔名與上一步的
※檔名第一個字元必須是字母，字母間不可留空格。
函數名稱相同即可。
8
• Step 6:編輯設定初始條件的M檔案
依Step 3再開啟一個新的M檔案，輸入紅色圈圈內的指
令。
下一頁會進行詳細說明
9
function setvalue
global a0 k timeup;
a0=input('請輸入反應物初始濃度:');
k=input('請輸入k值:');
timeup=input('請輸入時間區間上限:');
end
1. 第一行的部分不必像描述微分方程的M檔案那樣，只需要有函數
名稱即可，函數名稱可自訂，但不可與變數名稱相同，且第一個
字元必須為字母，字母間不可空格。
2. 第二行則是要將a0, k, timeup宣告為全域變數(global)，全域變數就
是要讓這些變數在各種檔案或模式中是互通共用的，因此這一步
驟是非常重要的。
3. 第三到五行則是利用input這個內建函數，讓使用者可以自行輸入
a0, k, timeup這些變數的值；紫色文字的部分會顯示給使用者看，
用來提醒使用者目前是要輸入哪個變數，紫色文字的部分可自行
更改。
4. 第四行的end跟第一行開頭的function是互相對應的， 代表該函數
運算的結束，不能漏掉。
10
• Step 7:將M檔案存檔
點選左上角 File → Save As，儲存檔案，檔名與上一步的
※檔名第一個字元必須是字母，字母間不可留空格。
函數名稱相同即可。
描述方程式及設定初始條件值的M檔案編輯完
成後，就可以進行MatLab內建的求解器函數進
行計算及作圖了。
11
Part II:命令視窗的計算與作圖
(以一級反應為例)
• Step 1:移至命令視窗
有>>出現的視窗即為命令視窗，命令視窗是供使用者輸入
函數、命令、運算式與程式語句的介面。
12
• Step 2:將變數宣告為全域變數
在命令視窗中輸入:
global a0 k timeup;
在命令視窗與M檔案中的變數原本是互相獨立
的，如果我們要使變數互通共用的話，我們必
須在M檔案以及命令視窗中都將a0,k,timeup宣
告為全域變數，若有一邊無宣告的話，則變數
依然會是獨立的，造成錯誤。
13
• Step 3:初始條件設定
執行Part I 所編輯好的setvalue.m，進行初始條件的設定。
直接在命令視窗中
輸入第二個M檔案
的檔名即可執行。
執行setvalue後，會請user
輸入參數值。這裡的值可
自行更改。
14
• Step 4:進行計算
這裡要空格
輸入以下指令:
[T,A]=ode45(@(t,a) firstorder(t,a,k),[0:0.1:timeup],a0);
下頁進行該指令詳細
說明
15
[T,A]=ode45(@(t,a) firstorder(t,a,k),[0:0.1:timeup],a0);
1. T和A這兩個變數是用來儲存計算後的結果，後面作圖也是利用這
兩個變數。
2. ode45是MatLab 內建的求解器，本次實驗使用的都是這個求解器。
3. @(t,a)是要告訴求解器:我們所要求的未知數只有 t 和 a。
4. 這裡就是先前M檔案中的函數名稱(輸入參數)，必須跟M檔案中一
致。
5. 這個部分是要設定計算的時間範圍，[開始時間:時間間隔:結束時
間] 。
6. 這個部分則是要設定初始條件，因為在這裡我們只求反應物的變
化，所以在這裡只要給定一個初始條件值 a0 即可，意思也就是
[A]0=a0；若要給定多個初始條件值時，這裡則必須用以下的表示
法:[a0 b0 c0]，每個初始條件以空格區隔開即可。
16
• Step 5:作圖
輸入plot(T,A)後，出現的Figure1就是一級反應time對濃度
所作的圖。
plot(X,Y)是作圖的函數，
第一個值為圖的X軸(Time)，
第二個值為圖的Y軸([A])。
這裡也可以直接輸入運算式
Ex: plot(T,A/a0)
產生的圖
X軸就是Time
Y軸就是[A]/[A]0
17
• Step 6:圖形設定
在Figure 1上選取Edit →Axes Properties，會出現紅色圈圈
的部分。
18
• Step 7:進行軸與圖形參數的微調，並且設定X軸Y軸以及
Title。
19
• Step 8:點選軸上的文字，可更改字型、大小、顏色等。
點選軸上文字或標題
可更改字體
大小及字型
20
• Step 9:將圖形存檔
圖形參數設定完之後，點選左上角File→Save As，儲存檔
案。
存檔類型若為*.fig，只有
MatLab能開啟。
若要一般圖檔則選取*.jpg。
21
• Step 10:抓取計算的Data
存完圖檔後，點選紅色圈圈處的Workspace，Workspace中
會存放著目前所有變數的值。
這裡會顯示變數的值或是矩陣大小，像A和T都是101X1
的行向量。
這裡的A和T矩陣大小一定要相同，若不同，則可能是
先前變數未清除完全，此時則必須先在命令視窗中輸入
clear all
(此指令會消除所有變數)
接著再重複執行Part II。
22
• Step 11:點選所要瀏覽的變數後，會出現下面的視窗，我們
便可以看到這個變數所有的值了。
23
• Step 12:選取所有數據，點滑鼠右鍵Copy，再開啟Excel，
將數據貼上，重複步驟將T和A的數據複製至Excel，也可
利用Excel作圖。
24
Part III:M檔案的撰寫
(以二級反應為例)
二級反應的微分式為:
d[A]/dt=-k*[A]^2
二級反應Part III的部分，與一級反應Part I的部分大同小異，僅有編輯
描寫微分方程的M檔案(第26頁)的地方有些微不同；若已熟悉Part I 的
同學，這裡可自行加快速度，若還不熟悉，可參考詳細步驟再練習一
次。
• Step 1:開啟Editor視窗 (可參考第4頁)
點選左上角File→New →Function，會出現一個Editor視窗。
25
• Step 2:編輯描寫微分方程的M檔案
在Editor視窗中輸入以下指令:
function da=secondorder(t,a,k)
da=zeros(1,1);
da(1)=-k*a(1)^2;
end
[這些指令下一頁會進行
詳細的說明]
26
function da=secondorder(t,a,k)
da=zeros(1,1);
da(1)=-k*a(1)^2;
end
二級反應的微分式為:
d[A]/dt=-k*[A]^2
1. 檔案的內容中，第一行必須遵循以下格式:
function d輸出參數=函數名稱(輸入參數)
function
da
= secondorder
(t,a,k)
這裡 a 是輸出參數。
而其中函數名稱可自訂 ，但不可與變數名稱相同，且第一個
字元必須為字母，字母間不可空格。
這裡 輸入參數共有三個，分別為 t, a & k 。
27
function da=secondorder(t,a,k)
da=zeros(1,1);
da(1)=-k*a(1)^2;
end
二級反應的微分式為:
d[A]/dt=-k*[A]^2
2. 第二行是將輸出參數“da “ 透過“zeros(1,1)” 將其設定為一個
1x1其值為零的矩陣。這個1x1的矩陣是要用來存放[A]值。所
以a(1)所代表的就是[A]。
3. 第三行則就用來描述二級反應的微分式， “ * ”號代表數學上的
乘號， “ ^ ”代表數學上的次方 。
4. 第四行的end跟第一行開頭的function是互相對應的，代表該函
數運算的結束，但不能漏掉。
※在MatLab中，於程式行末尾加上分號，會使該行程式執行結果不顯示。
28
• Step 3:將M檔案存檔
(可參考第8頁)
• Step 4:編輯設定初始條件的M檔案
(可參考第9、10
頁)
二級反應使用的條件設定M檔案與一級反應的相同，不熟
悉的同學可再練習一次。
29
Part IV:命令視窗的計算與作圖
(以二級反應為例)
二級反應Part IV的部分，與一級反應Part II的部分大同小異，僅有進
行計算時的指令(第31頁)有些微不同；若已熟悉Part II的同學，這裡
可自行加快速度，若還不熟悉，可參考詳細步驟再練習一次。
• Step 1:將變數宣告為全域變數 (可參考第13頁)
在命令視窗中輸入: global a0 k timeup;
• Step 2:初始條件設定
(可參考第14頁)
執行Part I 所編輯好的setvalue.m，進行初始條件的設定。
30
• Step 3:進行計算
這裡要空格
輸入以下指令:
[T,A]=ode45(@(t,a) secondorder(t,a,k),[0:0.1:timeup],a0);
下頁進行說明
31
[T,A]=ode45(@(t,a) secondorder(t,a,k),[0:0.1:timeup],a0);
1. T和A這兩個變數是用來儲存計算後的結果，後面作圖也是利用
這兩個變數。
2. ode45是MatLab內建的求解器，本次實驗使用的都是這個求解
器。
3. @(t,a)是要告訴求解器:我們所要求的未知數只有t和a。
4. 這裡就是先前M檔案中的函數名稱(輸入參數)，必須跟M檔案
中一致。
5. 這個部分是要設定計算的時間範圍，[開始時間:時間間隔:結束
時間] 。
6. 這個部分則是要設定初始條件，因為在這裡我們只求反應物的
變化，所以在這裡只要給定一個初始條件值a0即可，意思也就
是[A]0=a0；若要給定多個初始條件值時，這裡則必須用以下的
表示法:[a0 b0 c0]，每個初始條件以空格區隔開即可。
32
• Step 4:作圖 (可參考第17頁)
輸入plot(T,A)後，出現的Figure1就是二級反應time對 濃度
所作的圖。
• Step 5:圖形設定 (可參考第18頁)
在Figure 1上選取Edit →Axes Properties。
• Step 6:進行軸與圖形參數的微調，並且設定X軸Y軸以及
Title。
(可參考第19頁)
33
• Step 7:點選軸上的文字，可更改字型、大小、顏色等
(可參考第20頁)
• Step 8:將圖形存檔
(可參考第21頁)
圖形參數設定完之後，點選左上角File→Save As，儲存檔
案。
• Step 9~11:抓取計算的Data
(可參考第22~24頁)
存完圖檔後，將Workspace中T和A變數的值複製到Excel
並利用Excel作圖。
34
Part V:M檔案的撰寫
(以平行一級反應為例)
平行一級反應Part V的部分，與一級反應Part I的部分大同小異，僅有編
輯描寫微分方程的M檔案(第36頁)和編輯設定初始條件的M檔案(第
39頁)的地方有明顯不同；若已熟悉Part I 的同學，其他步驟可自行加
快速度，若還不熟悉，可參考詳細步驟再練習一次。
• Step 1:開啟Editor視窗 (可參考第4頁)
點選左上角File→New →Function，會出現一個Editor視
窗。
平行一級反應的微分式為:
AB
A C
反應速率常數= k1
反應速率常數= k2
d[A]/dt=-k1*[A]-k2*[A]
d[B]=k1*[A]
d[C]=k2*[A]
35
• Step 2:編輯描寫微分方程的M檔案
在Editor視窗中輸入以下指令:
function da=parallel(t,a,k1,k2)
da=zeros(3,1);
da(1)=-k1*a(1)-k2*a(1);
da(2)=k1*a(1);
da(3)=k2*a(1);
end
平行一級反應的微分式為:
d[A]/dt=-k1*[A]-k2*[A]
d[B]=k1*[A]
d[C]=k2*[A]
[這些指令下一頁會進行
詳細的說明]
36
function da=parallel(t,a,k1,k2)
da=zeros(3,1);
% a(1)=[A] a(2)=[B] a(3)=[C]
da(1)=-k1*a(1)-k2*a(1); ※在MatLab中，%是程式
da(2)=k1*a(1);
的註解符號，%後的指令
da(3)=k2*a(1);
文字皆不執行。
end
平行一級反應的微分式為:
d[A]/dt=-k1*[A]-k2*[A]
d[B]=k1*[A]
d[C]=k2*[A]
1. 檔案的內容中，第一行必須遵循以下格式:
function d輸出參數=函數名稱(輸入參數)
function
da
= parallel
(t,a,k)
這裡 a 是輸出參數。
而其中函數名稱可自訂 ，但不可與變數名稱相同，且第一個
字元必須為字母，字母間不可空格。
這裡 輸入參數共有三個，分別為 t, a & k 。
37
function da=parallel(t,a,k1,k2)
da=zeros(3,1);
% a(1)=[A] a(2)=[B] a(3)=[C]
da(1)=-k1*a(1)-k2*a(1); ※在MatLab中，%是程式
da(2)=k1*a(1);
的註解符號，%後的指令
da(3)=k2*a(1);
文字皆不執行。
end
平行一級反應的微分式為:
d[A]/dt=-k1*[A]-k2*[A]
d[B]=k1*[A]
d[C]=k2*[A]
2. 第二行的部分，因為計算後的值是利用行向量儲存的，
而這裡我們有[A], [B], [C]三個物種濃度，所以我們這裡
要將da設成(3,1)的矩陣。
3. 第三到第五行是用來描述平行一級反應的微分式， “ * ”號
代表數學上的乘號。
4. 第六行的end跟第一行開頭的function是互相對應的，代表該函
數運算的結束，但不能漏掉。
※在MatLab中，於程式行末尾加上分號，會使該行程式執行結果不顯示。
38
• Step 3:將M檔案存檔
(可參考第8頁)
點選左上角File→Save As，儲存檔案，檔名與上一步的函
數名稱相同即可。
• Step 4:編輯設定初始條件的M檔案
再開啟一個新的M檔案，輸入紅色圈圈內的指令。
下一頁會進行詳細說明
39
function setparallel
global a0 b0 c0 k1 k2 timeup;
a0=input(‘請輸入反應物的初始濃度:');
b0=input('請輸入產物B的初始濃度:');
c0=input('請輸入產物C的初始濃度:');
k1=input('請輸入k1值:');
k2=input('請輸入k2值:');
timeup=input('請輸入時間區間上限:');
end
1. 第一行的部分不必像描述微分方程的M檔案那樣，只需要有函數
名稱即可，函數名稱可自訂，但不可與變數名稱相同，且第一個
字元必須為字母，字母間不可能空格。
2. 第二行的部分，因為平行反應所要設定的條件較多，因此必須將
這些變數都宣告為全域變數(global)。
3. 第三到八行則是利用input這個內建函數，讓使用者可以自行輸入
這些變數的值；紫色文字的部分會顯示給使用者看，用來提醒使
用者目前是要輸入哪個變數，紫色文字的部分可自行更改。
4. 第九行的end跟第一行開頭的function是互相對應的，不能漏掉。
※在MatLab中，於程式行末尾加上分號，會使該行程式執行結果不顯示。
40
• Step 5:將M檔案存檔
(可參考第11頁)
點選左上角File→Save As，儲存檔案，檔名與上一步的函
數名稱相同即可。 ※檔名第一個字元必須是字母，字母間不可留空
格。
描述方程式及設定初始條件值的M檔案編輯完
成後，就可以進行MatLab內建的求解器函數進
行計算及作圖了。
41
Part VI:命令視窗的計算與作圖
(以平行一級反應為例)
平行一級反應Part VI的部分，與一級反應Part II的部分大同小異，僅有
進行計算時的指令(第46頁)有些微不同；若已熟悉Part I 的同學，其他
步驟可自行加快速度，若還不熟悉，可參考詳細步驟再練習一次。
• Step 1:將變數宣告為全域變數 (可參考第13頁)
在命令視窗中輸入:
global a0 b0 c0 k1 k2 timeup;
42
• Step 2:初始條件設定
執行Part V 所編輯好的setparallel.m，進行初始條件的設
定。
直接在命令視窗中
輸入M檔案的檔名
即可執行。
執行setvalue後，會請user
輸入參數值。這裡的值可
自行更改。
43
• Step 3:進行計算
輸入以下指令:
這裡要空格
[T,A]=ode45(@(t,a) parallel(t,a,k1,k2),[0:0.1:timeup],[a0 b0 c0]);
下頁進行說明
44
[T,A]=ode45(@(t,a) parallel(t,a,k1,k2),[0:0.1:timeup],[a0 b0 c0]);
1. T和A這兩個變數是用來儲存計算後的結果，後面作圖也是
利用這兩個變數。
2. ode45是MatLab內建的求解器，本次實驗使用的都是這個
求解器。
3. @(t,a)是要告訴求解器:我們所要求的未知數只有t和a。
4. 這裡就是先前M檔案中的函數名稱(輸入參數)，必須跟M
檔案中一致。
5. 這個部分是要設定計算的時間範圍，[開始時間:時間間隔:
結束時間] 。
6. 這個部分則是要設定初始條件，因為我們要求反應物A及
產物B、產物C三個物種的濃度變化，所以在這裡必須要給
定三個初始條件[A]0=a0 [B]0=b0 [C]0=c0。
45
• Step 4:作圖 (可參考第17頁)
輸入plot(T,A)後，出現的Figure1就是平行一級反應time對
濃度所作的圖。
• Step 5:圖形設定 (可參考第18頁)
在Figure 1上選取Edit →Axes Properties。
• Step 6:進行軸與圖形參數的微調，並且設定X軸Y軸以及
Title。
(可參考第19頁)
• Step 7:點選軸上的文字，可更改字型、大小、顏色等
(可參考第20頁)
46
• Step 8:點選Insert →Legend後，會出現紅色圈圈的部分。
47
• Step 9: 在
上按右鍵可更改設定，在圖上的字上點
選滑鼠兩下，可以更改文字。
完成時如右圖:
48
• Step 10:將圖形存檔
圖形參數設定完之後，點選左上角File→Save As，儲存檔
案。
存檔類型若為*.fig，只有
MatLab能開啟。
若要一般圖檔則選取*.jpg。
49
• Step 11:抓取計算的Data
存完圖檔後，點選紅色圈圈處的Workspace，Workspace中
會存放著目前所有變數的值。
這裡會顯示變數的值或是矩陣大小，這裡A的矩陣大小是
101X3的原因是:裡面一次存放了3個濃度([A] [B] [C])的
值，但是行的部分一樣是101，這代表[A] [B] [C]都分別有
101個數據，與T吻合。
50
• Step 12:點選所要瀏覽的變數後，會出現下面的視窗，我
們便可以看到這個變數所有的值了。
51
• Step 13:選取所有數據，點滑鼠右鍵Copy，再開啟Excel，
將數據貼上，重複步驟將T和A的數據複製至Excel，並利
用Excel作圖。
52