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(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)');