12/22/2006

12.8 非線性方程式解

12.8非線性方程式解


當常微分方程屬非線性時,就很難得到解析解,必須使用數值分析解。但這也常會發生一些垃圾答案的問題,必須就常識判斷其解之正當性。

例:有一大漏斗倉裝滿水,其外觀呈倒角錐形,半徑R高度為H,其體積為:

 V=πR²H/3

倉體之上方,有一進水口,其流入速率為Qin;其底部則另有孔口流出。截面積為A,其流量Q與倉中之水位高度h有下列關係:

Q = Cd A (2gh)½

方程式中,Cd為孔口係數,其值為0.6。根據流體質量不滅的定理,倉中體積的變化應等於流體進出之流量,即:

  dV/dt = Qin-Q

而由於r = (Rh/H),故:

 dV/dt=d/dt[πr²dh/3]=πR²h²/(3H²]dh/dt , 故

 [πR²h²/(3H²]dh/dt=Qin-Cd*A*(2gh)½  整理之,

dh/dt=[Qin-Cd*A*(2gh)½]/[πR²h²/(3H²]

整理上面之公式寫成程式,並加以執行。


function [h,t]=tri_tank(tf,h0,R,H,cd,A,Qin)
% Using Euler method to solve differential eqs.
% Water tank
% Inputs:
% h0:initial water level
% tf:time interval[t1 t2]
% R,H:radius & height of tank, m
% cd:coefficient of orifice
% A: Crosssection of orifice, m^2
% Qin: Inlet flow, cms
% Example: h=tri_tank([0 10],9,10,9,0.6,1,10)
aa=3*(H/R)^2/pi;g=9.8;
vv=@(t,h) aa*(Qin-cd*A*sqrt(2*g*h));
[t,h]=ode23(vv,tf,h0);
plot(t,h,'r-',t,h,'.');xlabel('t'),ylabel('y');

執行結果:

>>hold on; for i=[5 7 10 15],h=tri_tank([0 10],9,10,9,0.6,1,i);end


圖中所示,由上而下分別為Qin=15, 10, 7, 5 cms。前兩者之水位高度會累積比初設值高,後二者因進水量較低,故會消耗一部份水量,最後在4米高時達到平衡。

高階微分方程式


解二階以上之微分方程式必須採用特定的方法。例如考慮二階微方如下:

  5y"+7y'+4y = f(t)

先將其改寫為如下之型式以配合原來之一階微方:

  y"=(1/5)f(t)-(4/5)y-(7/5)y'


此時因為右邊仍有一階微分項,為令其符合一階之格式,可以定義另外新變數x1及x2:

x1=y
x2=y'


根據此兩定義項,可以再進行微分一次,得到下面之等式:


x1'=y'
x2'=y"


將此代入上項微方,得聯立方程式:


x1'=x2
x2'=y"=(1/5)f(t)-(4/5)x1-(7/5)x2


此型式稱為柯西式(Cauchy form) 。設 f(t) = sin(t),其區間為0<t<6 ,並設起始條件y(0)=3,y'(0)=9,輸入初值因此為[3, 9]。則可以利用ode23解法器求得相關值。由於變數x1及x2可以置於同一矩陣中,其函數因此可以呼叫如下:

fun=@(t,x) [x(2);(1/5)*sin(t)-4*x(1)-7*x(2)];
[t,x]=ode23(fun,[0 6],[3 9]);
plot(t,x(:,1),'r-',t,x(:,2),'bo',t,x(:,2))