MATLAB 程式設計 常微分方程 與聯立微分方程模擬 方煒 台大生機系

ODE 指令列表 MATLAB 用於求解常微分方程式的指令: 指令 方法 適用ODE類別 ode45 Explicit Runge-Kutta (4, 5) pair of Dormand-Prince Nonstiff ODE ode23 Explicit Runge-Kutta (2, 3) pair of Bogacki and Shampine ode113 Variable order Adams-Bashforth-Moulton PECE solver ode15s Numerical differentiation formulas (NDFS) Stiff ODE ode23s Modified Rosenbrock formula of order 2 ode23t Trapezoidal rule with a “free” interpolant ode23tb Implicit Runge-Kutta formula with a backward differentiation formula of order two

ODE 指令列表 指令主要可分兩大類 適用於 Nonstiff 系統 適用Stiff 系統 一般的常微分方程式都是 Nonstiff 系統 直接選用 ode45、ode23 或 ode113 來求解 適用Stiff 系統 速率（即微分值）差異相常大 使用一般的 ode45、ode23 或 ode113 來求解，可能會使得積分的步長（Step Sizes）變得很小，以便降低積分誤差至可容忍範圍以內，會導致計算時間過長 專門對付 Stiff 系統的指令，例如 ode15s、ode23s、ode23t 及 ode 23tb

ODE 指令基本用法 使用 ODE 指令時，必須先將要求解的 ODE 表示成一個函式 輸入為 t（時間）及 y（狀態變數，State Variables） 輸出則為 dy（狀態變數的微分值） ODE 函式的檔名為 odeFile.m，則呼叫 ODE 指令的格式如下： [t, y] = solver ('odeFile', [t0, t1], y0); [t0, t1] 是積分的時間區間 y0 代表起始條件（Initial Conditions） solver 是前表所列的各種 ODE 指令 t 是輸出的時間向量 y 則是對應的狀態變數向量

ODE 指令基本用法 以 van der Pol 微分方程為例，其方程式為： 化成標準格式 可用向量來表示成一般化的數學式 為一向量，代表狀態變數

ODE 指令基本用法 假設 =1， ODE 檔案（vdp1.m）可顯示如下： >> type vdp1.m function dy = vdp1(t, y) mu = 1; dy = [y(2); mu*(1-y(1)^2)*y(2)-y(1)]; 有了 ODE 檔案，即可選用前述之 ODE 指令來求解

範例-1 (I) 在 =1 時，van der Pol 方程式並非 Stiff 系統，所以使用ode45來畫出積分的結果 範例1：odeBasic01.m ode45('vdp1', [0 25], [3 3]');

範例-1 (II) 因為沒有輸出變數，所以上述程式執行結束後，MATLAB 只會畫出狀態變數對時間的圖形 [0, 25] 代表積分的時間區間，[3 3]’ 則代表起始條件（必須以行向量來表示） 因為沒有輸出變數，所以上述程式執行結束後，MATLAB 只會畫出狀態變數對時間的圖形

範例-2 (I) 要取得積分所得的狀態變數及對應的時間，可以加上輸出變數以取得這些資料 範例2：odeGetData01.m [t, y] = ode45('vdp1', [0 25], [3 3]'); plot(t, y(:,1), t, y(:,2), ':'); xlabel('Time t'); ylabel('Solution y(t) and y''(t)'); legend('y(t)', 'y''(t)');

範例-2 (II)

範例-3 (I) 也可以畫出 及 在 相位平面（Phase Plane ）的運動情況 範例3：odePhastPlot01.m [t, y] = ode45('vdp1', [0 25], [3 3]'); plot(y(:,1), y(:,2), '-o'); xlabel('y(t)'); ylabel('y''(t)');

範例-3 (II) 當 值越來越大時，van der Pol 方程式就漸漸變成一個 Stiff 的系統，此時若要解此方程式，就必須改用專門對付 Stiff 系統的指令

範例-4 (I) 將 值改成 1000， ODE 檔案改寫成（vdp2.m）： >> type vdp2.m function dy = vdp2(t, y) mu = 1000; dy = [y(2); mu*(1-y(1)^2)*y(2)-y(1)];

範例-4 (II) 選用專門對付 Stiff 系統的 ODE 指令，例如 ode15s，來求解此系統並作圖顯示 範例4：ode15s01.m [t, y]= ode15s('vdp2', [0 3000], [2 1]'); subplot(2,1,1); plot(t, y(:,1), '-o'); xlabel('Time t'); ylabel('y(t)'); subplot(2,1,2); plot(t, y(:,2), '-o'); xlabel('Time t'); ylabel('y''(t)'); % 注意單引號「'」的重覆使用

範例-4 (III) 由上圖可知， 的變化相當劇烈（超過 ），這就是 Stiff 系統的特色

範例-5 (I) 若要畫出二維平面相位圖，可以使用下列範例： 若要產生在某些特定時間點的狀態變數值，則呼叫 ODE 指令的格式可改成： 範例5：ode15s02.m 若要產生在某些特定時間點的狀態變數值，則呼叫 ODE 指令的格式可改成： [t, y] = solver('odeFile', [t0, t1, …, tn], y0); 其中 [t0, t1, …, tn] 即是特定時間點所形成的向量 [t, y]= ode15s('vdp2', [0 3000], [2 1]'); subplot(1,1,1); plot(y(:, 1), y(:, 2), '-o'); xlabel('y(t)'); ylabel('y''(t)')

範例-5 (II)

ODE 指令的選項 ODE 指令可以接受第四個輸入變數，代表積分過程用到的各種選項（Options），此種 ODE 指令的格式為： [t, y] = solver('odeFile', [t0, tn], y0, options); 其中 options 是由 odeset 指令來控制的結構變數 結構變數即包含了積分過程用到的各種選項 odeset 的一般格式如下： options = odeset('name1', value1, 'name2', value2, …) 其欄位 name1 的值為 value1，欄位 name2 的值為 value2，依此類推 未被設定的欄位，其欄位值即為預設值

ODE 指令的選項 也可以只改變一個現存的 options 結構變數中，某個欄位的值，其格式如下： newOptions = odeset(options, 'name', value); 若要讀出某個欄位的值，可用 odeget，其格式如下： value = odeget(otpions, 'name'); 其中 name 為欄位名稱，value 即為對應的欄位值 當使用 odeset 指令時，若無任何輸入變數，則 odeset 會顯示所有的欄位名稱及欄位值，並以大括號代表預設值

ODE 指令的選項 >> odeset AbsTol: [ positive scalar or vector {1e-6} ] RelTol: [ positive scalar {1e-3} ] NormControl: [ on | {off} ] NonNegative: [ vector of integers ] OutputFcn: [ function_handle ] OutputSel: [ vector of integers ] Refine: [ positive integer ] Stats: [ on | {off} ] InitialStep: [ positive scalar ] MaxStep: [ positive scalar ] BDF: [ on | {off} ] MaxOrder: [ 1 | 2 | 3 | 4 | {5} ] Jacobian: [ matrix | function_handle ] JPattern: [ sparse matrix ] Vectorized: [ on | {off} ] Mass: [ matrix | function_handle ] MStateDependence: [ none | {weak} | strong ] MvPattern: [ sparse matrix ] MassSingular: [ yes | no | {maybe} ] InitialSlope: [ vector ] Events: [ function_handle ]

由 odeset 產生的 ODE 選項 類別 欄位名稱 資料型態 預設值 說明 誤差容忍度之相關欄位 RelTol 正純量 相對誤差容忍度 AbsTo1 正純量或向量 絕對誤差容忍度 積分輸出之相關欄性 OutPutFcn 字串 ‘odeplot’ 輸出函式（若 ODE 指令無輸出變數，則在數值積分執行完畢後，MATLAB 會呼叫此輸出函式） OutputSel 索引向量 全部 ODE 指令之輸出變數的索引值，以決定那些輸出變數之元素將被送到輸出函式 Refine 正整數 1或 4 (for ode45) Refine = 2 可產生兩倍數量的輸出點，Refine = 3 可產生三倍數量的輸出點，依此類推。 Stats on 或 off off Stats = 'on' 產生計算過程的各種統計資料。

由 odeset 產生的 ODE 選項 Jacobian 矩陣之相關欄位 Jconstant on 或 off off ，則 Jacobian = 'on‘ Jpattern 是稀疏矩陣，則 JPattem = 'on' Vectorized 若 F(t, [y1, y2…..]) 傳回 [F(t,y1), F(t,y2)…..]，則 Vectorized = 'on' 積分步長（Step Sizes）之相關欄位 Max Step 正純量 ODE 指令之積分步長 的上限 Initial Step 起始步長的建議值 若F(t, y, ’Jacobian') 傳回 若 F([ ], [ ], ’JPattern’) 傳回 ，且

由 odeset 產生的 ODE 選項 質量矩陣之相關欄位 Mass none，M， M(t)，或 M(t, y) none MassSing ular yes，no 或 maybe 表明質量矩陣是否為 Singular 事件發生時間之相關欄位 Events on 或 off off 若 ODE 檔案並傳回事件函式及相關資訊，則 Events = 'on' Ode15s 之相關欄位 MaxOrder 1, 2, 3, 4 或 5 5 積分公式用到的最高次數 BDF 若使用 BDF（Backward Differentiation Formula）則 BDF = 'on'

常用欄位說明 f 在積分誤差容忍度方面，每一次積分所產生的局部誤差 e(i)，必須滿足下列方程式： max(RelTol* , AbsTol(i)) 其中i 代表第 i 個狀態變數 降低 RelTol 及 AbsTol 來求得更精確的積分結果

範例-1 (I) 範例6：odeRelTol01.m subplot(2,1,1); ode45('vdp1', [0 25], [3 3]'); title('RelTol=0.01'); options = odeset('RelTol', 0.00001); subplot(2,1,2); ode45('vdp1', [0 25], [3 3]', options); title('RelTol=0.0001');

範例-1 (II) 第一個圖所使用的相對誤差值是0.01（預設值），第二個圖所使用的相對誤差值是0.00001，因此我們得到較細密的點，但所花的計算時間也會比較長

積分輸出 積分輸出的相關處理方面 選用一個 OutputFcn 當 ODE 指令沒有輸出變數時，此輸出函式 OutputFcn 會被 MATLAB 呼叫 OutputFcn 的預設值是”odeplot”，其功能為畫出所有的狀態變數 其它可用的函式 odephas2：畫出 2-D 的相位平面（Phase Plane） odephas3：畫出 3-D 的相位平面 odeprint：隨時在指令視窗印出計算結果

積分輸出 以 Lorenz 渾沌方程式（Lorenz Chaotic Equation）為例 >> type lorenzOde.m function dy = lorenzOde(t, y) % LORENZODE: ODE file for Lorenz chaotic equation IGMA = 10.; RHO = 28; BETA = 8./3.; A = [ -BETA 0 y(2) 0 -SIGMA SIGMA -y(2) RHO -1 ]; dy = A*y;

範例-2 使用 ode45 對Lorenz 渾沌方程式進行數值積分 上圖中共有三條曲線，代表三個狀態變數隨時間變化 範例7：odeLorenz01.m 上圖中共有三條曲線，代表三個狀態變數隨時間變化 ode45('lorenzOde', [0 10], [20 5 -5]');

範例-3 上述範例畫三度空間之相位圖形 範例8：odeLorenz02.m 圖形中只出現一條曲線，此曲線代表以三個狀態變數為座標、以時間為參數的一條三度空間中的曲線 options = odeset('OutputFcn', 'odephas3'); % 使用 odephas3 進行繪圖 ode45('lorenzOde', [0 10], [20 5 -5]', options);

提示 要觀看 Lorenz 渾沌方程式隨時間而變的動畫，可在 MATLAB 指令視窗下直接輸入lorenz 指令

範例-4 (I) 假設 OutputFcn 設成“myfunc”： options = odeset('OutputFcn', 'myfunc') ODE 指令會呼叫 myfunc(tspan, y0, ‘init’) 讓 myfunc 進行各種初始化動作 積分步驟中，ODE 指令會持續呼叫status=myfunc(t, y) 若 status=1，則停止積分 積分結束時，ODE 指令會呼叫 myfunc([ ], [ ], ‘done’)，讓 myfunc 進行收尾動作 OutputSel 可指定要傳送到 OutputFcn 的狀態變數之元素

範例-4 (II) 只要傳送第一及第三個 Lorenz 渾沌方程式的狀態變數至 odeplot 範例9：odeOutputSelect01.m options = odeset('OutputSel', [1,3]) % 只畫出第一和第三個狀態變數 ode45('lorenzOde', [0 10], [20 5 -5]', options);

範例-5 (I) Refine 欄位可以使用內差法來增加輸出狀態變數的密度，以得到較平滑的輸出曲線 用 Refine 欄位使 ode23 的輸出點個數增為原先三倍： 範例10：odeRefine01.m subplot(2,1,1); ode23('vdp1', [0 25], [3 3]'); subplot(2,1,2); options = odeset('Refine', 3); ode23('vdp1', [0 25], [3 3]', options);

範例-5 (II)

範例-6 當 Stat=on 時，ODE 指令會在執行完畢後顯示計算過程的各種統計數字 範例11：odeShowStats01.m 71 successful steps 10 failed attempts 487 function evaluations [t, y] = ode45('vdp1', [0 25], [3 3]', odeset('Stat', 'on'));

範例-7 相同的統計數字，也可由 ODE 指令的第三個輸出變數傳回 範例12：odeShowStats02.m s = 71 10 487 [t, y, s] = ode45('vdp1', [0 25], [3 3]'); s

說明 MaxStep 及 InitialStep 欄位可用來調整最大積分步長及起始積分步長 一般而言，不必去調整這兩個數值，因為 ODE 指令本身就具有步長自動調適功能 若要產生更多輸出點，可直接調整 Refine 欄位值。調整 MaxStep 雖然可以達到同樣效果，但是計算時間可能會大幅增加 如果積分結果不甚準確，請勿先調降 MaxStep，您應先調降 RelTol 及 AbsTol。調降 MaxStep 是最後的步驟

進階用法 ODE 檔案的進階用法，使 ODE 指令能夠迅速且準確地算出積分結果 可將 tspan（積分時間範圍）、y0（起始值）及 options（ODE參數）置於 ODE 檔案中，這些變數必須能由 ODE 檔案傳回，其格式為： [tspan, y0, options] = odeFile([], [], 'init') 假設 odeFile 即是我們的 ODE 檔案且 odeFile 滿足上述要求，則可以直接呼叫 ODE 指令如下： [t, y] = solver('odeFile') 其中 solver 為前述的任一個 ODE 指令，它可由 odeFile 直接得到 tspan、y0 及 options 等積分所需的資訊

範例-1 (I) 以前述的 van der Pol 為例，若要能夠傳回 tspan、y0 及 options，vdp1.m 須改寫如下（vdp3.m）： >> type vdp3.m function [output1, output2, output3] = vdp3(t, y, flag) if strcmp(flag, '') mu = 1; output1 = [y(2); mu*(1-y(1)^2)*y(2)-y(1)]; % dy/dt elseif strcmp(flag, 'init'), output1 = [0; 25]; % Time span output2 = [3; 3]; % Initial conditions output3 = odeset; % ODE options end

範例-1 (II) 範例13：odeAdvanced01.m ode45('vdp3')

範例-2 (I) van der Pol 的微分方程式有一個參數 ，希望從外面傳入此參數的值(vdp4.m ) >> type vdp4.m function [output1, output2, output3] = vdp4(t, y, flag, mu) if nargin < 4 | isempty(mu), mu = 1; end if strcmp(flag, '') output1 = [y(2); mu*(1-y(1)^2)*y(2)-y(1)]; % dy/dt elseif strcmp(flag, 'init'), output1 = [0; 25]; % Time span output2 = [3; 3]; % Initial conditions output3 = odeset; % ODE options

範例-2 (II) 就變成一個選擇性（Optional）的參數，其預設值為 1 將 的值從 MATLAB 傳入，並畫出不同 值下的 van der Pol 方程式的狀態變數： 範例14：odeAdvanced02.m subplot(2,1,1); ode45('vdp4', [], [], [], 1); % mu=1 subplot(2,1,2); ode45('vdp4', [], [], [], 3); % mu=3

範例-2 (III) 用到了許多空矩陣，這些空矩陣代表「取用預設值」，因此 ode45 會直接從 vdp4.m 取用時間區間及變數起始值 的值分別是 1 及 3 用到了許多空矩陣，這些空矩陣代表「取用預設值」，因此 ode45 會直接從 vdp4.m 取用時間區間及變數起始值 也可以傳入二個或更多的參數，MATLAB 及 ODE 指令對於可傳入的參數個數並無設限

進階用法列表 為解決其它較複雜的 ODEs 及 DAEs（Differential Algebra Equations），ODE 檔案亦可在不同的旗號（Flag）下傳回不同的資訊，列表如下： 旗號 傳回值 （空字串） dy(=F(t,y)) init tspan, yo 及 options jacobian Jacobian 矩陣 J(t,y) = 2F/2Y jpattern Jacobian sparsity pattern 之矩陣 mass 解 M(t,y)y' = F(t,y) 所須的 質量矩陣 M events 定義事件發生點的各種資訊

Write the ODE file (I) for y’’’ + y’’ + y’ =0 using ‘odefuc1.m’ as filename. Step 1: define y1=y, y2 = y’, y3=y’’ Step 2: expressions: y’1=y2, y’2=y3, y’3=y’’’=-y’’-y’=-y3-y2

Write the ODE file (II) % odefuc1.m function ydot= odefuc1(t,y); %ydot=[y(2);y(3);-y(3)-y(2)]; T1=y(2); T2=y(3); T3=-y(3)-y(2); ydot=[T1;T2;T3];

Please check ode-sup.pdf for more details. 補充 1 % odesup1.m clear; clc y0=[10;1;0]; tspan=[0 20]; [t,y]=ode45('odefuc1',tspan,y0); subplot(3,1,1); plot(t,y(:,1)); xlabel('t'); ylabel('dy/dt'); subplot(3,1,2); plot(t,y(:,2)); xlabel('t'); ylabel('d^2y/dt^2'); subplot(3,1,3); plot(t,y(:,3)); xlabel('t'); ylabel('d^3y/dt^3'); Please check ode-sup.pdf for more details.

補充 2 % odesup2.m clear; clc; tic y0=[10;1;0]; tspan=[0 20]; [t,y,s]=ode45('odefuc1',tspan,y0); toc plot(t,y(:,1)); xlabel('t'); ylabel('dy/dt'); s

補充 3 % odesup3.m clear; clc; y0=[10;1;0]; tspan=[0 20]; options=odeset(‘stats’,’on’); [t,y]=ode45('odefuc1',tspan,y0,options); plot(t,y(:,1)); xlabel('t'); ylabel('dy/dt');

補充 4 % odesup4.m clear; clc; tic y0=[10;1;0]; tspan=[0 20]; options=odeset(‘stats’,’on’,’RelTol’,1e-2); [t,y]=ode23s('odefuc1',tspan,y0,options); toc plot(t,y(:,1)); xlabel('t'); ylabel('dy/dt');

補充 5 % odesup5.m clear; clc y0=[10;1;0]; tspan=[0 20]; options=odeset('OutputFcn','odephas3','OutputSel',[1 2 3]); [t,y]=ode45('odefuc1',tspan,y0,options); xlabel('y(1)');ylabel('y(2)'); zlabel('y(3)');

補充 6 % odesup6.m clear; clc y0=[10;1;0]; tspan=[0 20]; subplot(2,2,1); options=odeset('OutputFcn','odephas2','OutputSel',[1 2]); [t,y]=ode45('odefuc1',tspan,y0,options); xlabel('y(1)');ylabel('y(2)'); subplot(2,2,2); options=odeset('OutputFcn','odephas2','OutputSel',[2 3]); xlabel('y(2)');ylabel('y(3)'); subplot(2,2,3); options=odeset('OutputFcn','odephas2','OutputSel',[1 3]); xlabel('y(1)');ylabel('y(3)'); subplot(2,2,4); options=odeset('OutputFcn','odephas2','OutputSel',[3 1]); xlabel('y(3)');ylabel('y(1)');

補充 7 % odesup7.m clear; clc y0=[10;1;0]; tspan=[0 20]; subplot(2,1,1); options=odeset('OutputFcn','odephas3','OutputSel',[1 2 3]); [t,y]=ode45('odefuc1',tspan,y0,options); xlabel('y(1)');ylabel('y(2)'); zlabel('y(3)'); subplot(2,1,2); options=odeset('OutputFcn','odephas3','OutputSel',[1 2 3],'refine',20); xlabel('y(1)');ylabel('y(2)'); zlabel('y(3)');

聯立微分方程模擬

CUC01.m % Temperature regime in the soil layer CUC01.m % Boundary condition is surface temperature. % Function required: soil01.m % function cuc01(trial) global ks cs if nargin==0 % If no argument trial=1; % use 1 as the default. end if trial>4 | trial <1 % If argument_value >4 or <1 trial =1 % use 1 as the default switch trial case 1 ks=5.5; cs=2000; case 2 ks=11; cs=2000; % ks doubled case 3 ks=5.5; cs=4000; % cs doubled case 4 ks=11; cs=4000; % ks, cs both doubled % ks: Soil thermal conductivity (kJ/m/C) and ks/3.6 (W/m/C) % cs: Heat capacity of soil (kJ/m3/C)

補充 8 (III) % tstart = 0; tfinal = 48; y0 = [10; 10; 10; 10; 10] ; % A 5 by 1 matrix for initial conditions. [t,y] = ode23t('soil01',[tstart tfinal],y0); % calling function ode23t with constants and eqs. in 'soil01.m' % with simulated time from tstart to tfinal % with initial conditions in matrix y0 and % with calculated answer in matrix y. plot(t,y(:,1),'b^-',t,y(:,2),'gV-',t,y(:,3),'r+-',t,y(:,4),'c*-',t,y(:,5),'ko-'); axis([-inf,inf,5,15]); grid on; xlabel('time elapsed, hr'); ylabel('Soil temperature, ^oC'); tit=['Given conditions: ks=',num2str(ks),' and cs=', num2str(cs)]; title(tit); legend('T1','T2','T3','T4','T5',2); fprintf('\n Running Case %1.0f: given ks=%3.0f and cs=%4.0f. \n\n', trial,ks, cs); disp(' Please check Figure_Window for simulated results.'); disp(' Totally, 5 curves showing T1, T2, T3, T4 and T5 versus t.');

Soil01.m % soil01.m % function dy = soil01(t,y) global ks cs z=0.1; % Depth of each soil layer (m) T0=10; % Average outside temperature (C) TU=5; % Amplitude, range of temperature variation (C) TBL=10; % Boundary soil temperature (C) TF=T0+TU*sin(2.*pi/24.*(t-8.)); % t is time in hours % TF: Soil temperature of surface layer (C) % Maximum temperature occurs at 2 o'clock in the afternoon T1=y(1);T2=y(2);T3=y(3);T4=y(4);T5=y(5); DIFF_T1=2*(TF-T1)+(T2-T1); DIFF_T2=(T1-T2)+(T3-T2); DIFF_T3=(T2-T3)+(T4-T3); DIFF_T4=(T3-T4)+(T5-T4); DIFF_T5=(T4-T5)+(TBL-T5)*2; val=ks/z/z/cs; dy = [DIFF_T1; DIFF_T2; DIFF_T3; DIFF_T4; DIFF_T5]*val; % Format of dy needs to be consistent with y0 in cuc01.m that is a 5 by 1 matrix % end of function

The End