12.8 非線性方程式解

12.8非線性方程式解

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

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

　V=πR²H/3

倉體之上方，有一進水口，其流入速率為Ｑin；其底部則另有孔口流出。截面積為Ａ，其流量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。前兩者之水位高度會累積比初設值高，後二者因進水量較低，故會消耗一部份水量，最後在４米高時達到平衡。

高階微分方程式

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

　　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))