1 / 52

使用 MatLab 解微分方程

使用 MatLab 解微分方程. 使用 MatLab 解微分方程式的基本流程與相關基本知識. 首先撰寫描述該微分方程式的 M 檔案。 M 檔案就像是一個小程式，裡面可以是欲執行的命令、函數、運算式或程式語句。 將該 M 檔案與初始條件值，載入 ode 求解函數中，透過求解器求出微分方程的解。 最後經過 ode 求解函數運算後，將傳回指定的時間區間內的解，以一個 行向量 儲存。 行向量是一個 列數恆為 1 的矩陣，矩陣在 MatLab 中表示為 ( 行 , 列 ) ， 例如 : (1,1) 、 (2,1) …

quinn-hicks
Download Presentation

使用 MatLab 解微分方程

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. 使用 MatLab解微分方程

  2. 使用 MatLab 解微分方程式的基本流程與相關基本知識 • 首先撰寫描述該微分方程式的M檔案。M檔案就像是一個小程式，裡面可以是欲執行的命令、函數、運算式或程式語句。 • 將該M檔案與初始條件值，載入ode求解函數中，透過求解器求出微分方程的解。 • 最後經過ode求解函數運算後，將傳回指定的時間區間內的解，以一個行向量儲存。 • 行向量是一個列數恆為1的矩陣，矩陣在MatLab中表示為(行,列) ， 例如:(1,1)、(2,1) … (1,1)代表”1 x 1”的矩陣，可以儲存一個變數值； (2,1)代表”2 x 1”的矩陣，可以儲存兩個變數值；以此類推。 5. 因此要解微分方程時，我們必須先撰寫好兩個M檔案，一個是用來描述該微分方程式，一個則是用來進行初始條件值的設定。

  3. Part I: M檔案的撰寫(以一級反應為例) 一級反應的微分式為: d[A]/dt = -k*[A] • Step 1:開啟MatLab程式 點擊桌面上的MatLab圖示，執行MatLab。 (MatLab剛開啟時較慢，請耐心等候。) • Step 2:待紅色圈圈處的>>及Ready出現，即為開啟完成。

  4. Step 3:開啟 Editor視窗 點選左上角 File→New →Function，會出現一個 Editor視窗。

  5. Step 4:編輯描寫微分方程的M檔案 在Editor視窗中輸入以下指令: function da=firstorder(t,a,k) da=zeros(1,1); da(1)=-k*a(1); end [這些指令下一頁會進行詳細的說明]

  6. function da=firstorder(t,a,k) da=zeros(1,1); da(1)=-k*a(1); end 一級反應的微分式為: d[A]/dt=-k*[A] • 檔案的內容中，第一行必須遵循以下格式: function d輸出參數=函數名稱(輸入參數) function da = firstorder (t,a,k) 這裡 a 是輸出參數。 而其中函數名稱可自訂，但不可與變數名稱相同，且第一個 字元必須為字母，字母間不可空格。 這裡 輸入參數共有三個，分別為 t, a & k 。

  7. function da=firstorder(t,a,k) da=zeros(1,1); da(1)=-k*a(1); end 一級反應的微分式為: d[A]/dt=-k*[A] • 第二行是將輸出參數“da“透過“zeros(1,1)”將其設定為一個1x1其值為零的矩陣。這個1x1的矩陣是要用來存放[A]值。所以a(1)所代表的就是[A]。 • 第三行則就用來描述一級反應的微分式， “*”號代表數學上的乘號。 • 第四行的end跟第一行開頭的function是互相對應的，代表該函數運算的結束，但不能漏掉。 ※在MatLab中，於程式行末尾加上分號，會使該行程式執行結果不顯示。

  8. ※檔名第一個字元必須是字母，字母間不可留空格。※檔名第一個字元必須是字母，字母間不可留空格。 • Step 5:將M檔案存檔 點選左上角 File→Save As，儲存檔案，檔名與上一步的函數名稱相同即可。

  9. Step 6:編輯設定初始條件的M檔案 依Step 3再開啟一個新的M檔案，輸入紅色圈圈內的指令。 下一頁會進行詳細說明

  10. functionsetvalue global a0 k timeup; a0=input('請輸入反應物初始濃度:'); k=input('請輸入k值:'); timeup=input('請輸入時間區間上限:'); end 第一行的部分不必像描述微分方程的M檔案那樣，只需要有函數名稱即可，函數名稱可自訂，但不可與變數名稱相同，且第一個字元必須為字母，字母間不可空格。 第二行則是要將a0, k, timeup宣告為全域變數(global)，全域變數就是要讓這些變數在各種檔案或模式中是互通共用的，因此這一步驟是非常重要的。 第三到五行則是利用input這個內建函數，讓使用者可以自行輸入a0, k, timeup這些變數的值；紫色文字的部分會顯示給使用者看，用來提醒使用者目前是要輸入哪個變數，紫色文字的部分可自行更改。 第四行的end跟第一行開頭的function是互相對應的， 代表該函數運算的結束，不能漏掉。

  11. ※檔名第一個字元必須是字母，字母間不可留空格。※檔名第一個字元必須是字母，字母間不可留空格。 • Step 7:將M檔案存檔 點選左上角 File → Save As，儲存檔案，檔名與上一步的函數名稱相同即可。 描述方程式及設定初始條件值的M檔案編輯完成後，就可以進行MatLab內建的求解器函數進行計算及作圖了。

  12. Part II:命令視窗的計算與作圖(以一級反應為例) • Step 1:移至命令視窗 有>>出現的視窗即為命令視窗，命令視窗是供使用者輸入函數、命令、運算式與程式語句的介面。

  13. Step 2:將變數宣告為全域變數 在命令視窗中輸入:global a0 k timeup; 在命令視窗與M檔案中的變數原本是互相獨立的，如果我們要使變數互通共用的話，我們必須在M檔案以及命令視窗中都將a0,k,timeup宣告為全域變數，若有一邊無宣告的話，則變數依然會是獨立的，造成錯誤。

  14. Step 3:初始條件設定 執行Part I所編輯好的setvalue.m，進行初始條件的設定。 直接在命令視窗中輸入第二個M檔案的檔名即可執行。 執行setvalue後，會請user輸入參數值。這裡的值可自行更改。

  15. 這裡要空格 • Step 4:進行計算 輸入以下指令: [T,A]=ode45(@(t,a) firstorder(t,a,k),[0:0.1:timeup],a0); 下頁進行該指令詳細說明

  16. [T,A]=ode45(@(t,a) firstorder(t,a,k),[0:0.1:timeup],a0); T和A這兩個變數是用來儲存計算後的結果，後面作圖也是利用這兩個變數。 ode45是MatLab 內建的求解器，本次實驗使用的都是這個求解器。 @(t,a)是要告訴求解器:我們所要求的未知數只有 t和 a。 這裡就是先前M檔案中的函數名稱(輸入參數)，必須跟M檔案中一致。 這個部分是要設定計算的時間範圍，[開始時間:時間間隔:結束時間]。 這個部分則是要設定初始條件，因為在這裡我們只求反應物的變化，所以在這裡只要給定一個初始條件值 a0即可，意思也就是[A]0=a0；若要給定多個初始條件值時，這裡則必須用以下的表示法:[a0 b0 c0]，每個初始條件以空格區隔開即可。

  17. 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

  18. Step 6:圖形設定 在Figure 1上選取Edit →Axes Properties，會出現紅色圈圈的部分。

  19. Step 7:進行軸與圖形參數的微調，並且設定X軸Y軸以及Title。

  20. Step 8:點選軸上的文字，可更改字型、大小、顏色等。 點選軸上文字或標題 可更改字體大小及字型

  21. Step 9:將圖形存檔 圖形參數設定完之後，點選左上角File→Save As，儲存檔案。 存檔類型若為*.fig，只有MatLab能開啟。 若要一般圖檔則選取*.jpg。

  22. Step 10:抓取計算的Data 存完圖檔後，點選紅色圈圈處的Workspace，Workspace中會存放著目前所有變數的值。 這裡會顯示變數的值或是矩陣大小，像A和T都是101X1的行向量。 這裡的A和T矩陣大小一定要相同，若不同，則可能是先前變數未清除完全，此時則必須先在命令視窗中輸入 clear all(此指令會消除所有變數) 接著再重複執行Part II。

  23. Step 11:點選所要瀏覽的變數後，會出現下面的視窗，我們便可以看到這個變數所有的值了。

  24. Step 12:選取所有數據，點滑鼠右鍵Copy，再開啟Excel，將數據貼上，重複步驟將T和A的數據複製至Excel，也可利用Excel作圖。

  25. Part III:M檔案的撰寫(以二級反應為例) 二級反應的微分式為: d[A]/dt=-k*[A]^2 • Step 1:開啟Editor視窗(可參考第4頁) 點選左上角File→New →Function，會出現一個Editor視窗。 二級反應Part III的部分，與一級反應Part I的部分大同小異，僅有編輯描寫微分方程的M檔案(第26頁)的地方有些微不同；若已熟悉Part I 的同學，這裡可自行加快速度，若還不熟悉，可參考詳細步驟再練習一次。

  26. Step 2:編輯描寫微分方程的M檔案 在Editor視窗中輸入以下指令: function da=secondorder(t,a,k) da=zeros(1,1); da(1)=-k*a(1)^2; end [這些指令下一頁會進行詳細的說明]

  27. function da=secondorder(t,a,k) da=zeros(1,1); da(1)=-k*a(1)^2; end 二級反應的微分式為: d[A]/dt=-k*[A]^2 • 檔案的內容中，第一行必須遵循以下格式: function d輸出參數=函數名稱(輸入參數) function da = secondorder (t,a,k) 這裡 a 是輸出參數。 而其中函數名稱可自訂，但不可與變數名稱相同，且第一個 字元必須為字母，字母間不可空格。 這裡 輸入參數共有三個，分別為 t, a & k 。

  28. function da=secondorder(t,a,k) da=zeros(1,1); da(1)=-k*a(1)^2; end 二級反應的微分式為: d[A]/dt=-k*[A]^2 • 第二行是將輸出參數“da“透過“zeros(1,1)”將其設定為一個1x1其值為零的矩陣。這個1x1的矩陣是要用來存放[A]值。所以a(1)所代表的就是[A]。 • 第三行則就用來描述二級反應的微分式， “*”號代表數學上的乘號， “^”代表數學上的次方。 • 第四行的end跟第一行開頭的function是互相對應的，代表該函數運算的結束，但不能漏掉。 ※在MatLab中，於程式行末尾加上分號，會使該行程式執行結果不顯示。

  29. Step 3:將M檔案存檔(可參考第8頁) • Step 4:編輯設定初始條件的M檔案(可參考第9、10頁) 二級反應使用的條件設定M檔案與一級反應的相同，不熟悉的同學可再練習一次。

  30. Part IV:命令視窗的計算與作圖(以二級反應為例) 二級反應Part IV的部分，與一級反應Part II的部分大同小異，僅有進行計算時的指令(第31頁)有些微不同；若已熟悉Part II的同學，這裡可自行加快速度，若還不熟悉，可參考詳細步驟再練習一次。 • Step 1:將變數宣告為全域變數(可參考第13頁) 在命令視窗中輸入:global a0 k timeup; • Step 2:初始條件設定(可參考第14頁) 執行Part I所編輯好的setvalue.m，進行初始條件的設定。

  31. 這裡要空格 • Step 3:進行計算 輸入以下指令: [T,A]=ode45(@(t,a) secondorder(t,a,k),[0:0.1:timeup],a0); 下頁進行說明

  32. [T,A]=ode45(@(t,a) secondorder(t,a,k),[0:0.1:timeup],a0); T和A這兩個變數是用來儲存計算後的結果，後面作圖也是利用這兩個變數。 ode45是MatLab內建的求解器，本次實驗使用的都是這個求解器。 @(t,a)是要告訴求解器:我們所要求的未知數只有t和a。 這裡就是先前M檔案中的函數名稱(輸入參數)，必須跟M檔案中一致。 這個部分是要設定計算的時間範圍，[開始時間:時間間隔:結束時間]。 這個部分則是要設定初始條件，因為在這裡我們只求反應物的變化，所以在這裡只要給定一個初始條件值a0即可，意思也就是[A]0=a0；若要給定多個初始條件值時，這裡則必須用以下的表示法:[a0 b0 c0]，每個初始條件以空格區隔開即可。

  33. Step 4:作圖 (可參考第17頁) 輸入plot(T,A)後，出現的Figure1就是二級反應time對濃度 所作的圖。 • Step 5:圖形設定(可參考第18頁) 在Figure 1上選取Edit →Axes Properties。 • Step 6:進行軸與圖形參數的微調，並且設定X軸Y軸以及Title。(可參考第19頁)

  34. Step 7:點選軸上的文字，可更改字型、大小、顏色等 (可參考第20頁) • Step 8:將圖形存檔(可參考第21頁) 圖形參數設定完之後，點選左上角File→Save As，儲存檔案。 • Step 9~11:抓取計算的Data (可參考第22~24頁) 存完圖檔後，將Workspace中T和A變數的值複製到Excel 並利用Excel作圖。

  35. Part V:M檔案的撰寫(以平行一級反應為例) 平行一級反應Part V的部分，與一級反應Part I的部分大同小異，僅有編輯描寫微分方程的M檔案(第36頁)和編輯設定初始條件的M檔案(第39頁)的地方有明顯不同；若已熟悉Part I 的同學，其他步驟可自行加快速度，若還不熟悉，可參考詳細步驟再練習一次。 • Step 1:開啟Editor視窗(可參考第4頁) 點選左上角File→New →Function，會出現一個Editor視窗。 平行一級反應的微分式為: d[A]/dt=-k1*[A]-k2*[A] d[B]=k1*[A] d[C]=k2*[A] A  B 反應速率常數= k1 A C 反應速率常數= k2

  36. 平行一級反應的微分式為: d[A]/dt=-k1*[A]-k2*[A] d[B]=k1*[A] d[C]=k2*[A] • 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 [這些指令下一頁會進行詳細的說明]

  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); 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] ※在MatLab中，%是程式的註解符號，%後的指令文字皆不執行。 • 檔案的內容中，第一行必須遵循以下格式: function d輸出參數=函數名稱(輸入參數) function da = parallel (t,a,k) 這裡 a 是輸出參數。 而其中函數名稱可自訂，但不可與變數名稱相同，且第一個 字元必須為字母，字母間不可空格。 這裡 輸入參數共有三個，分別為 t, a & k 。

  38. 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); 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] ※在MatLab中，%是程式的註解符號，%後的指令文字皆不執行。 • 第二行的部分，因為計算後的值是利用行向量儲存的，而這裡我們有[A],[B], [C]三個物種濃度，所以我們這裡要將da設成(3,1)的矩陣。 • 第三到第五行是用來描述平行一級反應的微分式， “*”號代表數學上的乘號。 • 第六行的end跟第一行開頭的function是互相對應的，代表該函數運算的結束，但不能漏掉。 ※在MatLab中，於程式行末尾加上分號，會使該行程式執行結果不顯示。

  39. Step 3:將M檔案存檔(可參考第8頁) 點選左上角File→Save As，儲存檔案，檔名與上一步的函數名稱相同即可。 • Step 4:編輯設定初始條件的M檔案 再開啟一個新的M檔案，輸入紅色圈圈內的指令。 下一頁會進行詳細說明

  40. functionsetparallel global a0 b0 c0 k1 k2 timeup; a0=input(‘請輸入反應物的初始濃度:'); b0=input('請輸入產物B的初始濃度:'); c0=input('請輸入產物C的初始濃度:'); k1=input('請輸入k1值:'); k2=input('請輸入k2值:'); timeup=input('請輸入時間區間上限:'); end 第一行的部分不必像描述微分方程的M檔案那樣，只需要有函數名稱即可，函數名稱可自訂，但不可與變數名稱相同，且第一個字元必須為字母，字母間不可能空格。 第二行的部分，因為平行反應所要設定的條件較多，因此必須將這些變數都宣告為全域變數(global)。 第三到八行則是利用input這個內建函數，讓使用者可以自行輸入這些變數的值；紫色文字的部分會顯示給使用者看，用來提醒使用者目前是要輸入哪個變數，紫色文字的部分可自行更改。 第九行的end跟第一行開頭的function是互相對應的，不能漏掉。 ※在MatLab中，於程式行末尾加上分號，會使該行程式執行結果不顯示。

  41. Step 5:將M檔案存檔(可參考第11頁) 點選左上角File→Save As，儲存檔案，檔名與上一步的函數名稱相同即可。 ※檔名第一個字元必須是字母，字母間不可留空格。 描述方程式及設定初始條件值的M檔案編輯完成後，就可以進行MatLab內建的求解器函數進行計算及作圖了。

  42. Part VI:命令視窗的計算與作圖(以平行一級反應為例) • Step 1:將變數宣告為全域變數(可參考第13頁) 在命令視窗中輸入:globala0b0 c0 k1 k2 timeup; 平行一級反應Part VI的部分，與一級反應Part II的部分大同小異，僅有進行計算時的指令(第46頁)有些微不同；若已熟悉Part I 的同學，其他步驟可自行加快速度，若還不熟悉，可參考詳細步驟再練習一次。

  43. Step 2:初始條件設定 執行Part V所編輯好的setparallel.m，進行初始條件的設定。 直接在命令視窗中輸入M檔案的檔名即可執行。 執行setvalue後，會請user輸入參數值。這裡的值可自行更改。

  44. 這裡要空格 • Step 3:進行計算 輸入以下指令: [T,A]=ode45(@(t,a) parallel(t,a,k1,k2),[0:0.1:timeup],[a0 b0 c0]); 下頁進行說明

  45. [T,A]=ode45(@(t,a) parallel(t,a,k1,k2),[0:0.1:timeup],[a0 b0 c0]); T和A這兩個變數是用來儲存計算後的結果，後面作圖也是利用這兩個變數。 ode45是MatLab內建的求解器，本次實驗使用的都是這個求解器。 @(t,a)是要告訴求解器:我們所要求的未知數只有t和a。 這裡就是先前M檔案中的函數名稱(輸入參數)，必須跟M檔案中一致。 這個部分是要設定計算的時間範圍，[開始時間:時間間隔:結束時間]。 這個部分則是要設定初始條件，因為我們要求反應物A及產物B、產物C三個物種的濃度變化，所以在這裡必須要給定三個初始條件[A]0=a0 [B]0=b0 [C]0=c0。

  46. 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頁)

  47. Step 8:點選Insert →Legend後，會出現紅色圈圈的部分。

  48. 在 • Step 9:上按右鍵可更改設定，在圖上的字上點 選滑鼠兩下，可以更改文字。 完成時如右圖:

  49. Step 10:將圖形存檔 圖形參數設定完之後，點選左上角File→Save As，儲存檔案。 存檔類型若為*.fig，只有MatLab能開啟。 若要一般圖檔則選取*.jpg。

  50. Step 11:抓取計算的Data 存完圖檔後，點選紅色圈圈處的Workspace，Workspace中會存放著目前所有變數的值。 這裡會顯示變數的值或是矩陣大小，這裡A的矩陣大小是101X3的原因是:裡面一次存放了3個濃度([A][B][C])的值，但是行的部分一樣是101，這代表[A][B][C]都分別有101個數據，與T吻合。

More Related