12/22/2006

12.7 實例-剛性與非剛性問題

(編譯及圖示摘自mathworks.com)

簡單非剛性問題

MATLAB提供一個名為rigidode之程式是處理一剛性體體在無外在受力之情況下,所產生之速度與位移變化。這是一個由柯洛氏(Krogh)針對非剛性解法器提出的標準測試方法。其解析之答案屬於 Jacobian楕圓函數,可以由MATLAB中得到。其所用的區間約為1.5週期。其求解之聯立方程式如下:



為解此聯立微分方程,必須先建立包含這些方程式之函數,設為rigid,其中y為代表聯立方程式右邊之函數關係,為含有三元素之行向量,其內容如下:

function dydt = rigid(t,y)
dydt =[y(2) * y(3);-y(1) * y(3);-0.51 * y(1) * y(2)];


上述函數亦可利用隱函數表示,如下:

rigid=@(t,y) [y(2) * y(3);-y(1) * y(3);-0.51 * y(1) * y(2)];


解法器可用ode45,其呼叫格式為:ode45(rigid,,),其中y0為初值,即為[0 1 1];tspan經歷時間,此處設為[0 12]。為方便起見,整個問題之定義及解題以M-檔存放。其微分方程則用次函數f表示。由於程式呼叫ode45 指令函數,沒有輸出項,故其輸出屬內定輸出函數odeplot,利用該函數直接將結果繪圖。


function rigidodex(tspan,y0)
% Solution of Euler's rigid body Equation
% revised from rigidode.m
% tspan: Interval, [t1, t2]
% y0: Initial value
%Example: rigidodex([0 12],[0 1 1])
y0 =y0(:);
rigid=@(t,y) [y(2)*y(3);-y(1)*y(3);-0.51*y(1)*y(2)];
options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);
[T,Y]=ode45(rigid,tspan,y0); % solve the problem using ODE45
plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')


執行上述程式之結果如下圖:



剛性問題


凡德波(van der Pol)方程式屬一般剛性問題,其剛性用µ表示(本例之預設值設為µ=1000)。其微分聯立方程式如下:



在第二式中,當值增加,問題的本質會更顯出其剛性,此時擺動週期也隨之增大。當其值為1000時,擺動進入遲滯狀態,其剛性更強。在其限制的週期內有一部份解變化趨緩,然後反向迅速變化,呈現交替狀態。初期值則接近於變化緩慢之區域,以測試步距之大小。

為解上述聯立方程式,可先建立對應之函數vdp1000:

function dy = vdp1000(t,y)
dy = zeros(2,1); % a column vector
dy(1) = y(2);
dy(2) = 1000*(1 - y(1)^2)*y(2) - y(1);


對於剛性問題所用之ODE解法器係以採用近似Jacobian矩陣為主。為此通常可用odeset('Jacobian',@J)進行設定。其中J為次函數J(t,y,mu),屬於上項係數之偏微分矩陣,可求得點(t,y)附近具µ值之Jacobian 矩陣df/dx。

本例中,初值為[2 0],時間區間為[0 3000],使用ode15s解法器求解。

[T,Y] = ode15s(@vdp1000,[0 3000],[2 0]);

完整的程式如下:

function vdpodex(MU,y0)
% Solution of van der Pol Equation
% Example: vdpodex(500,[2 0])
if nargin==0,MU=1000;y0=[2 0];end
if nargin==1,;y0=[2 0];end
y0 =y0(:);tspan=[0; max(20,3*MU)];
J=@(t,y,mu) [0 1;-2*mu*y(1)*y(2)-1 mu*(1-y(1)^2) ];
f=@(t,y,mu) [y(2) ; mu*(1-y(1)^2)*y(2)-y(1) ];
options = odeset('Jacobian',J);
figure;
[t,y]=ode15s(f,tspan,y0,options,MU);
plot(t,y(:,1),'-o');


執行結果如下圖:



µ=1,並將繪圖指令改為:

plot(t,y(:,1),'-',t,y(:,2),'-.')

則上述程式執行後情況略有不同,其結果如下:

>>vdpodex(1,[2 1]) %設mu=1



具時間變數之微分方程


本例旨在解具有時間變數項之微分方程式。以下面之 ODE為例,其時間因變參數僅經由一系列之資料以兩向量型式決定:

y'(t) + f(t)y(t) = g(t)

式中之初值條件為 y(0) = 0,其中函數 f(t)定義為 n-x-1 向量 tf 與 f;而函數 g(t)則經過 m-x-1 向量 tg 與 g定義。

首先,與時間變化有關之參數 f(t)與 g(t)可以表示如下:

ft = linspace(0,5,25); % Generate t for f
f = ft.^2 - ft - 3; % Generate f(t)
gt = linspace(1,6,25); % Generate t for g
g = 3*sin(gt-0.25); % Generate g(t)

寫出 M-檔案函數程式就上述之指定值進行內插,以求得特定時間點之時變對應值。

function dydt = myode(t,y,ft,f,gt,g)
f = interp1(ft,f,t); % Interpolate the data set (ft,f) at time t
g = interp1(gt,g,t); % Interpolate the data set (gt,g) at time t
dydt = -f.*y + g; % Evalute ODE at time t

其次,利用ode45解法器中之參數呼叫上述之導數函數 myode.m :

Tspan = [1 5]; % Solve from t=1 to t=5
IC = 1; % y(t=0) = 1
[T Y] = ode45(@(t,y) myode(t,y,ft,f,gt,g),TSPAN,IC); % Solve ODE

繪出 y(t)解與時間之函數關係:

plot(T, Y);
title('Plot of y as a function of time');
xlabel('Time'); ylabel('Y(t)');