12/22/2006

12.6 一階微分方程之解法

12.6一階微分方程之解法


連續性函數之微分可以利用多種指令為之,而基本之細段微分法則是最早使用的概念,亦可利用MATLAB運算指令直接求解。

尤拉法(Euler method)


尤拉法為數值積分法最基本之一種。它是利用時間細段微分,利用初值及一階微分方程進行求解。以下列方程式為例:

y'=dy/dt = f(t,y)


式中t為時間,y為因變數。y', dy/dt 均為其一階微分的表示式。m為t之函數。根據微分之定義,可以改寫如下:

y'=dy/dt = limΔ-0 [y(t+Δt)-y(t)]/dt = f(t, y)

展開,可以求得經過Δt之時間後,y(t+Δt)之結果為:

y(t+Δt)=y(t)+Δtf(t,y)

若設Δt為一均勻之區段時間,或稱為步進值,則可以給予y(t)一個初值,然後利用增加步進值之個數逐步累積至最終的區間,最後得到y之答案值。其求解之型式如下:

y(tk+1) = y(tk) + Δtf(tk,y(tk))

故只要給定時間範圍及f(t0,y0)等初值,即可利用迴圈的方式得到最終答案。其中,Δt設定值之大小則會影響答案之準確度,Δt值愈大,準確度會差,Δt值愈小則運算時間加長。

例:設f(t,y)=-10y,初值t=0, y0=2。且y'=-10y 之真實解為y(t)=2e-10t。設Δt=0.02,以此利用MATLAB求解。

function y=euler(y0,time,delta)
% Using Euler method to solve differential eqs.
% Inputs:
% y0:initial value
% time:time limit
% delta:step size
% Example: y=euler(2,0.5,0.02)
r=-10;k=0;
t=0:delta:time;y=zeros(size(t));
y(1)=y0;
for i=2:length(t),y(i)=y(i-1)+r*y(i-1)*delta;end
y_true=2*exp(-10*t);
plot(t,y,'o',t,y_true),xlabel('t'),ylabel('y');

執行結果:



>> y=euler(2,0.5,0.02);




同樣的問題亦可使用ode45指令計算,則先設定其函數為fun:

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



利用MATLAB之ODE45之指令所得之結果比前面所用的方法準確度高。以下我們將介紹數種微分方程之解法及相關指令。