12/22/2006

12.7初值型微分方程

12.7初值型微分方程


常微分方程式(Ordinary differential equations, ODE)為一般工數中求解的問題,包含有自變數與應變數,可以應用在多工程多方面動態之解析。若自變數增多,則成偏微分方程式,簡稱PDE(Partial differential equation)。目前MATLAB提供之微分方程解法器有數種,其性能及應用範圍略有不同。其中較常用的有ode23及ode45,這兩個均是利用Runge-Kutta法進行求解。其餘如ode113、ode15s、ode15i、ode23s、ode23t、ode23tb等等,其功能與範圍在運用上均有些差異,有些是根據不同方法得到的求解過程,常應用在控制方面。

常微分初值型解法器之相關函數均是利用數值積分方法作為求解過程。最初由設定之起始時間與條件,依設定之時間驅間逐步計算其解。若結果能滿足解法器設定之容許標準,就算成功;若無法達到,則必須再降低區間重新嚐試。常微分解法器之輸出入格式如下:

 [t,Y] = solver(odefun,tspan,y0)
 [t,Y] = solver(odefun,tspan,y0,options)
 [t,Y,TE,YE,IE] = solver(odefun,tspan,y0,options)
 sol = solver(odefun,[t0 tf],y0...)

上述指令中之solver,所指的就是ode45, ode23, ode113, ode15s, ode23s, ode23t, or ode23tb等等。其他相關參數介紹如下:

odefun:


用以計算微分方右邊之函數。所有解法器處理之系統方程式均是y'=dydt=f(t,y)之型式,有些是包括質量矩之問題,如M(t,y)y'=f(t,y)之函數。這些都與時間t有關,t為常數,而dydt及y均為行向量。ode45s解法器僅能處理質量為常數之問題;而ode15s與ode23t則能處理具奇數型質量矩陣之方程式,或稱為Differential-algebraic equations(DAEs)。

tspan 積分區間之向量[t0, tf]。解法器設定初值在tspan(1),然後由tspan(1)積分至tspan(end)。為得到特定時間(升冪或降冪)之對應解,可以使用tspan=[t0,t1, . . .,tf]。若tspan=[t0, tf],則會回應每個積分之步驟值。若tspan多於二個元素,則會回應時間對應值。時間必須依順序安排,或升冪或降冪。由於內部運算之安排,使用tspan向量之大小並不影響運算時間,但會佔據較多的記憶體。

y0:初值向量。

options:改變積分特性之相關參數。這些參數也可用odeset函數取得或設定。

輸出值包括t,Y,前者為時間點之行向量;後者為解答陣列,每一行之解y均與時間t相對應。每組(t,Y)之產生稱為事件函數。每次均會檢查是否函數等於零。並決定是否在零時終止運算。這可以在函數中之特性上設定。例如以events 或@events產生一函數。

[value, isterminal,direction]=events(t,y)

其中,value(i)為函數之值,isterminal(i)=1時運算在等於零時停止,=0時繼續;direction(i)=0時所有零時均需計算(預設值), +1在事件函數增加時等於零, -1在事件函數減少時等於零等狀況。

此外,TE, YE, IE則分別為事件發生之時間,事件發生時之答案及事件函數消失時之指標i。

有關初值型常微分方程之指令及其功能說明如下:


指令名稱 、說明及算法:
  • ode45 非剛性中階解法器 Runge-Kutta
  • ode23 非剛性低階解法器 Runge-Kutta
  • ode113 非剛性變階解法器 Adams
  • ode15s 剛性變階解法器及 DAE NDFs (BDFs)
  • ode23s 剛性低階解法器及 DAE Rosenbrock
  • ode23t 中剛性梯形法則解法獸器及DAE Trapezoidal rule
  • ode23tb 剛性低階解法器 TR-BDF2
茲就這些指令名稱說明如下:

ode45


此指令以 Runge-Kutta (4,5) 演算法為基礎,運算上是屬單步驟計算y(tn),故僅需時間點前之解y(tn-1)。 ode45函數指令可作為第一階段評估的數據,然後再思考是否採取下一步驟。

ode23


與ode45函數同以 Runge-Kutta (2,3)方程式為基礎,並以Bogacki 與Shampine配對。在溫和剛性的情況,容許度較粗略時可能比ode45 更有效率。此指令也是單步驟解題。

ode15s


以數值微分方程(NDFs).為基礎之可變順序指令。可選擇後向微分方程BDFs, (或稱為Gear 法)。ode15s 是一個多步驟指令。若執行前認為問題特性屬於剛性,亦或使用過ode45 指令但無法達成目的時,可以試用f ode15s指令。

ode23s


以Rosenbrock修正方程式二階為基礎。. 由於它是單步驟解題指令,故在粗略容許範圍下可能比ode15s 更有效率。它可能解一些ode15s效率低或無法解的剛性問題。

ode23t


採用自由間極的梯形法。在問題屬於中等剛性或不需涵括數值阻尼的特性時可以使用此指令解題。

ode23tb


執行 TR-BDF2—初期利用梯形法涵蓋Runge-Kutta 方程式而第二階段採用後向微分法二階。 如同ode23s一樣,在粗略容許度下此法比ode15s有效率。

ODE 處理


積分指令之特性、參數可以藉ODE函數之處理而改變指令之內涵,並因而可以影響解題的結果:

函數名稱及說明
  • odeset 設定或改變ODE指令之輸入選項
  • odeget 取得odeset產生的選項特性

ODE 指令之輸出函數:


若輸出函數己經設定,則解題指令會在完成積分步驟後,呼叫這些特定函數。此時你可以使用odeset指令指定其中之一個樣本函數,作為OutputFcn 的特性,或你也可以加以修正,使其成為自已擁有的函數。

  函數名稱及說明
  • odeplot   時序圖
  • odephas2  二維平面圖
  • odephas3  三維相面圖
  • odeprint  印出至指令窗

例:設一微分方程式為y'=sin(t),其初值y(0)=0,區間為 ,設其步進為週期(2π)之1/13,或者Δt=2π/13。其實際積分值為1-cos(t)。以尤拉法寫出之程式如下:

function y=euler2(y0,time,n)
% Using Euler method to solve differential eqs.
% y'=A*sin(t)
% Inputs:
% y0:initial value
% time:time limit,in the form of [t1 t2]*pi
% n:in the form of [1/n]*t2*pi, for step size
% Example: y=euler2(0,[0 2],13)
k=0;delta=time(2)/n;tt=time*pi;
t=tt(1):delta:tt(2);y=zeros(size(t));
y(1)=y0;
for i=2:length(t),
y(i)=y(i-1)+sin(t(i-1))*delta;
end
y_true=1-cos(t);
plot(t,y,'o',t,y_true),xlabel('t'),ylabel('y');


執行結果:

>>y=euler2(0,[0 4],13);



若改用ode45指令,其程式如下:

fun=@(t,y) sin(t);
[t,y]=ode45(fun,[0 4*pi],0);
plot(t,y,'ro',t,1-cos(t)),xlabel('t'),ylabel('y');




顯然利用ode45解法器比利用迴圈法為佳,當然迴圈法之Δt設定越小值,其準確度會趨近於後者之狀況。在上式之fun函數數,注意即使沒有y變數參與其中,也要將其列入,否則無法執行。
例:試解下面之微分方程式


 10y' + y = te-2t y(0)=2; 0≦t≦2


並以其結果 y(t)=[732e-0.1t-19e-2t-10e-2t]/361比較之。
解:

fun=@(t,y) 0.1*(t*exp(-2*t)-y);
[t,y]=ode23(fun,[0 2],2);
y2=(732*exp(-0.1*t)-19*t.*exp(-2*t)-10*exp(-2*t))/361;
plot(t,y,'ro',t,y2),xlabel('t'),ylabel('y');



電路分析


例:一電路由一電阻器R與一電容器C串聯接於一電流v上,試求通過電容器兩端之端電壓y。
解:根據克希荷夫電壓定律總電壓為電流端電壓之和,故可建立如下之微分式:

RC[dy/dt]+y=v(t)
y(0)=2V RC=0.1s


由於開始時並無電源,故v=0。設時間常數RC=0.1s,經整理,其結果如次:

dy/dt =[v(t)-y]/0.1 = -10y


執行時,使用ode23解法器,一樣可獲得正確的結果。

fun=@(t,y) -10*y;
[t,y]=ode23(fun,[0, 0.4],2);
plot(t,y,'ro',t,2*exp(-10*t)),xlabel('t'),ylabel('y');



就所用點數而言,ode23解法器所用的點數較少。

由以上之範例可知,線性微分方程式均可使用數值法進行解析,雖然也有通用解可用,但大部份仍以數值法較為簡便。但若階數大於二時,使用通式會相當複雜,只有數值法可以解決,而且得到的結果可以立即用圖去表示,因而可以作比較。

若上例中之外在電壓有值,其型式為:

v(t) = 10et/Isin(2πt/p)

同樣,設 RC=0.1s,則前例仍可維持如下之型式:

dy/dt = 10[v(t)-y]


解法:

function y=euler3(y0,tf,tau,P)
% Using Euler method to solve differential eqs.
% RC(dy/dt)+y=v(t)
% Inputs:
% y0:initial value
% tf:time interval[t1 t2]
% tau, P:constants
% Example: y=euler3(0,[0 2],0.3,2)
vv=@(t) 10*exp(-t/tau).*sin(2*pi*t/P);
fun=@(t,y) 10*(vv(t)-y);
[t,y]=ode23(fun,tf,y0);
subplot(2,1,1)
plot(t,vv(t),'k-');xlabel('t'),ylabel('y');
subplot(2,1,2)
plot(t,y,'r-',t,y,'.');xlabel('t'),ylabel('y');

執行結果:

>>y=euler3(0,[0 2],0.3,2); %tau=0.3s, P=2s




>>y=euler3(0,[0 0.2],0.05,0.03); %tau=0.05s, P=0.03s



>>y=euler3(0,[0 0.8],0.3,0.03); %tau=0.3s, P=0.03s