| Mathematics | |
例題：ODE 初期値問題ソルバの適用
この節では、MATLAB の中で解くことのできる種々の種類の例題を示しています。
rigidode)vdpode)fem1ode)brussode)ballode)orbitode)hb1dae)例題：単純なノンスティッフ問題
rigidodeは、ノンスティッフソルバに対して、Krogh
[6]により提唱された標準テスト問題の解を示します。
ユーザの便利さのため、すべての問題は、一つの
M ファイルで定義され、解かれます。微分方程式は、サブ関数fとしてコード化されています。ode45
ソルバは、出力引数なしでコールすると、デフォルトの出力関数odeplotが使われて、解成分をプロットします。この例題を実行するには、コマンドラインで、単にrigidodeをタイプしてください。
この例題をMATLAB Help
ブラウザから実行するには、デモ名をクリックしてください。他の方法としては、コマンドラインで、rigidode と入力してください。
function rigidode
%RIGIDODE 外力のない剛体のEuler方程式
tspan = [0 12];
y0 = [0; 1; 1];
% ode45を使って、問題を解く
ode45(@f,tspan,y0);
% -------------------------------------------------------------
function dydt = f(t,y)
dydt = [ y(2)*y(3)
-y(1)*y(3)
-0.51*y(1)*y(2) ];
例題：スティッフな問題 (van der Pol 方程式)
vdpodeは、例題：van
der Pol 方程式、 µ = 1000 (Stiff)の中で記述されている van der Pol方程式を示します。微分方程式
この問題は、
が増加すると、よりスティッフになり、振動の周期が長くなります。
が1000
を超える場合、方程式は、過減衰になり、問題は強いスティッフ性を示します。リミットサイクルは、解の要素がゆっくり変化する部分と問題が非常にスティッフな部分を含んでいます。すなわち、スティッフでない部分と急激に変化する部分が交互に生じます。
デフォルトでは、スティッフな問題に対して使用するODE
suite の中のソルバは、Jacobian
行列を数値的に近似します。しかし、ここの例題では、サブ関数J(t,y,mu)
を与えて、
= muに対して、(t,y)で、Jacobian行列
を解析的に計算します。解析的なJacobianの使用は、積分の信頼性と効率を向上します。
MATLAB Helpブラウザから、この例題を実行するには、デモ名をクリックする、コマンドラインで、vdpode
と入力してください。コマンドラインで、vdpodeへの引数として、
の値を指定することができます。デフォルトは、
= 1000です。
を含んでいます。 function vdpode(MU)
%VDPODE パラメータ化した van der Pol 方程式(大きな MUに対して、スティッフ)
if nargin < 1
MU = 1000; % デフォルト
end
tspan = [0; max(20,3*MU)]; % 数周期
y0 = [2; 0];
options = odeset('Jacobian',@J);
[t,y] = ode15s(@f,tspan,y0,options,MU);
plot(t,y(:,1));
title(['Solution of van der Pol Equation, \mu = ' num2str(MU)]);
xlabel('time t');
ylabel('solution y_1');
axis([tspan(1) tspan(end) -2.5 2.5]);
---------------------------------------------------------------
function dydt = f(t,y,mu)
dydt = [ y(2)
mu*(1-y(1)^2)*y(2)-y(1) ];
---------------------------------------------------------------
function dfdy = J(t,y,mu)
dfdy = [ 0 1
-2*mu*y(1)*y(2)-1 mu*(1-y(1)^2) ];
例題：有限要素離散化
fem1odeは、偏微分方程式の有限要素離散からの結果の
ODE 問題の解を示します。fem1ode(N)の中のNの値は、離散化をコントロールするもので、求まる結果の方程式は、N個の方程式から構成されます。デフォルトは、N
= 19 です。
この例題は、質量行列を含んでいます。ODEs のシステムは、偏微分方程式のライン解法から求まります。
初期条件
と境界条件
を使います。整数
が選択され、
は
と定義されます。そして、偏微分方程式の解は、k = 0,
1, ..., N+1 に対して、
で、つぎの式
によって、近似されます。ここで、
は、
で1、すべての他の
で0である区分的線形関数です。Galerkin
離散化は、つぎの ODEs システムを導きます。
初期値
は、偏微分方程式に対して、初期条件から求まります。問題は、時間間隔
で解かれます。
と
は、つぎの関係があります。
options = odeset('Mass',@mass,'MStateDep','none','Jacobian',J)
は、問題が
の型をしていることを示しています。サブ関数 A(t,N)は、時変質量行列
を計算し、R
は、定数のJacobian になります。
MATLAB Help
ブラウザからこの例題を実行するには、デモ名をクリックするか、コマンドラインで、fem1ode
と入力してください。コマンドラインで、fem1odeの引数として、
の値を設定することができます。デフォルトは、
= 19です。
function fem1ode(N)
%FEM1ODE は、時変質量行列をもつスティッフな問題
if nargin < 1
N = 19;
end
h = pi/(N+1);
y0 = sin(h*(1:N)');
tspan = [0; pi];
% ヤコビアンは定数
e = repmat(1/h,N,1); % e=[(1/h) ... (1/h)];
d = repmat(-2/h,N,1); % d=[(-2/h) ... (-2/h)];
J = spdiags([e d e], -1:1, N, N);
options = odeset('Mass',@mass,'MStateDependence','none', ...
'Jacobian',J);
[t,y] = ode15s(@f,tspan,y0,options,N);
surf((1:N)/(N+1),t,y);
set(gca,'ZLim',[0 1]);
view(142.5,30);
title(['Finite element problem with time-dependent mass ' ...
'matrix, solved by ODE15S']);
xlabel('space ( x/\pi )');
ylabel('time');
zlabel('solution');
%--------------------------------------------------------------
-
function out = f(t,y,N)
h = pi/(N+1);
e = repmat(1/h,N,1); % e=[(1/h) ... (1/h)];
d = repmat(-2/h,N,1); % d=[(-2/h) ... (-2/h)];
J = spdiags([e d e], -1:1, N, N);
out = J*y;
%--------------------------------------------------------------
-
function M = mass(t,N)
h = pi/(N+1);
e = repmat(exp(-t)*h/6,N,1); % e(i)=exp(-t)*h/6
e4 = repmat(4*exp(-t)*h/6,N,1);
M = spdiags([e e4 e], -1:1, N, N);
例題 ：大規模、スティッフなスパース行列
brussodeは、(ポテンシャル的に)大規模なスティッフスパース問題の解を示します。問題は、古典的な"Brusselator"システム[2] で、化学反応の拡散モデルです。
と共に解かれます。システム内には、
個の方程式が存在し、方程式が、
であるように順序付けられている場合、Jacobian
は、制約幅 5をもつ帯域になります。
brussode(N) をコール中、パラメータN 
2がグリッド点の数を設定するのに使います。ここで、N
は、
に対応します。結果のシステムは、2N個の方程式から構成されます。デフォルトでは、N
は20です。問題は、Nが大きくなると、スティッフ性が強くなり、Jacobian
のスパース性が増大することです。
サブ関数F(t,y,N) は、Brusselator
問題用の微係数ベクトルを出力します。サブ関数 jpattern(N)は、Jacobian
の中の非ゼロの位置を示す1と0のスパース行列を出力します。デフォルトでは、スティッフな
ODE ソルバは、フル行列として、数値的な Jacobian
を生成します。プロパティJPattern は、ソルバが Jacobian
をスパース行列として、数値的に生成するために使われるスパースパターンをもったソルバを与えるために使われます。スパースパターンを与えると、Jacobian
を作成するために必要な関数計算の回数が著しく減り、シミュレーションが速くなります。Brusselator
問題に対して、スパース性が与えられなければ、関数の 2N
回の計算が、2N 行 2N 列の Jacobian
行列を計算するために必要となります。スパース性が、与えられた場合、N
の値に関わらず、四つの計算のみが必要です。
MATLAB Help
ブラウザからこの例題を実行するには、デモ名をクリックするか、コマンドラインで、brussode
と入力してください。コマンドラインで、brussodeへの引数として、
の値を設定することができます。デフォルトは、
= 20です。
= 1/50 とで
function brussode(N)
% BRUSSODE 化学反応をモデル化したスティッフな問題
if nargin<1
N = 20;
end
tspan = [0; 10];
y0 = [1+sin((2*pi/(N+1))*(1:N)); repmat(3,1,N)];
options = odeset('Vectorized','on','JPattern',jpattern(N));
[t,y] = ode15s(@f,tspan,y0,options,N);
u = y(:,1:2:end);
x = (1:N)/(N+1);
surf(x,t,u);
view(-40,30);
xlabel('space');
ylabel('time');
zlabel('solution u');
title(['The Brusselator for N = ' num2str(N)]);
% --------------------------------------------------------------
function dydt = f(t,y,N)
c = 0.02 * (N+1)^2;
dydt = zeros(2*N,size(y,2)); % dy/dtを前もって設定
% グリッドの片端で、関数の2成分を計算
i = 1;
dydt(i,:) = 1 + y(i+1,:).*y(i,:).^2 - 4*y(i,:) + ...
c*(1-2*y(i,:)+y(i+2,:));
dydt(i+1,:) = 3*y(i,:) - y(i+1,:).*y(i,:).^2 + ...
c*(3-2*y(i+1,:)+y(i+3,:));
% すべての内部グリッド点の関数の2成分を計算
i = 3:2:2*N-3;
dydt(i,:) = 1 + y(i+1,:).*y(i,:).^2 - 4*y(i,:) + ...
c*(y(i-2,:)-2*y(i,:)+y(i+2,:));
dydt(i+1,:) = 3*y(i,:) - y(i+1,:).*y(i,:).^2 + ...
c*(y(i-1,:)-2*y(i+1,:)+y(i+3,:));
% グリッドの他の片端で、関数の2成分を計算
i = 2*N-1;
dydt(i,:) = 1 + y(i+1,:).*y(i,:).^2 - 4*y(i,:) + ...
c*(y(i-2,:)-2*y(i,:)+1);
dydt(i+1,:) = 3*y(i,:) - y(i+1,:).*y(i,:).^2 + ...
c*(y(i-1,:)-2*y(i+1,:)+3);
% --------------------------------------------------------------
function S = jpattern(N)
B = ones(2*N,5);
B(2:2:2*N,2) = zeros(N,1);
B(1:2:2*N-1,4) = zeros(N,1);
S = spdiags(B,-2:2,2*N,2*N);
例題：単純イベント位置検出問題
ballode
は、バウンドするボールの運動をモデル化します。この例題は、ODE
ソルバがイベントの生じる位置を検出する機能を示すものです。
この例題で、イベント関数は、サブ関数eventsでコード化されています。
[value,isterminal,direction] = events(t,y)
-1 |
負の方向のみのゼロ交点を検出 |
0 |
すべてのゼロ交点の検出 |
1 |
正の方向のみのゼロ交点を検出 |
value, isterminal, directionの長さは、event
関数の数と同じになります。各ベクトルの i
番目の要素は、i
番目の関数に対応します。より進んだイベントの発生位置を検出する例題は、例題 :アドバンスなイベント位置検出問題を参照してください。
ballodeの中で、Eventsプロパティを@events
に設定することは、ボールが落ちる途中(direction = -1)で、ボールが地面にヒットした(高さ
y(1) は 0になる)とき、シミュレーションを停止(isterminal
= 1)
することになります。それで、例題では、バウンドするボールに対応する初期条件を付けて、積分を再スタートします。
この例題を、MATLAB Help
ブラウザから実行するには、デモ名をクリックするか、コマンドラインで、ballodeと入力してください。
function ballode
% BALLODE は、バウンドするボールのデモを実行します。
tstart = 0;
tfinal = 30;
y0 = [0; 20];
refine = 4;
options = odeset('Events',@events,'OutputFcn', @odeplot,...
'OutputSel',1,'Refine',refine);
set(gca,'xlim',[0 30],'ylim',[0 25]);
box on
hold on;
tout = tstart;
yout = y0.';
teout = [];
yeout = [];
ieout = [];
for i = 1:10
% 最初のターミナルイベントまで解析
[t,y,te,ye,ie] = ode23(@f,[tstart tfinal],y0,options);
if ~ishold
hold on
end
% 出力をまとめる
nt = length(t);
tout = [tout; t(2:nt)];
yout = [yout; y(2:nt,:)];
teout = [teout; te]; % tstartでのイベントをレポートしない。
yeout = [yeout; ye];
ieout = [ieout; ie];
ud = get(gcf,'UserData');
if ud.stop
break;
end
% .9 の減衰をもつ新しい初期条件を設定
y0(1) = 0;
y0(2) = -.9*y(nt,2);
% 最初の timestep の良い推定は、最後の timestep の長さになり、それを
% 使うおとで、計算は速くなります。
options = odeset(options,'InitialStep',t(nt)-t(nt-refine),...
'MaxStep',t(nt)-t(1));
tstart = t(nt);
end
plot(teout,yeout(:,1),'ro')
xlabel('time');
ylabel('height');
title('Ball trajectory and the events');
hold off
odeplot([],[],'done');
% --------------------------------------------------------------
function dydt = f(t,y)
dydt = [y(2); -9.8];
% --------------------------------------------------------------
function [value,isterminal,direction] = events(t,y)
% 高さが、減じる方向でゼロを通る時刻を求め、シミュレーションを
% 停止する。
value = y(1); % 高さ = 0 を検出
isterminal = 1; % シミュレーションを停止
direction = -1; % 負の方向のみ
例題：アドバンスなイベント位置検出問題
orbitodeは、ノンスティッフな問題に対して使用するソルバに対する標準テストの解を示します。これは、月の周りを回り、地球に戻ってくる宇宙船の軌跡を示すものです(Shampine
and Gordon [6], p.246)。
最初の2つの解要素は、無限小の質量の座標系で、他のものに対してプロットしたものが、他の2つの物体の周りを回る軌跡になります。初期条件は、軌跡が周期性をもつように選択されます。
の値は、宇宙船が月の周りを回り、地球に戻るように設定します。中程度の厳密さをもったトレランスが、軌跡の質的に高い挙動を再生するために必要となります。適切な値は、RelTol が1e-5で、AbsTolが1e-4です。
eventsサブ関数は、スタート点から最大の距離をもつ点の位置を示すイベント関数と宇宙船が、スタート点に戻る時刻を含むものです。積分に使われるステップサイズが大きくて、イベントの位置を決定できない場合でさえ、イベントの正確な位置を割り出すことに注意してください。この例題で、ゼロクロッシングの方向を決定する機能は、クリティカルなものです。初期の点と最大距離の点として出力される点は共に、同じイベント関数値をもち、クロッシングの方向は、それらを区別するために使われます。
MATLAB Help
ブラウザからこの例題を実行するには、デモ名をクリックするか、コマンドラインで、orbitode
と入力してください。この例題は、出力関数 odephase2 を使って、2次元の位相平面プロットを作成し、積分の進行状況を見ることができます。
function orbitode
% ORBITODE は、制約付き3体問題です。
tspan = [0 7];
y0 = [1.2; 0; 0; -1.04935750983031990726];
options = odeset('RelTol',1e-5,'AbsTol',1e-4,...
'OutputFcn',@odephas2,'Events',@events);
[t,y,te,ye,ie] = ode45(@f,tspan,y0,options);
plot(y(:,1),y(:,2),ye(:,1),ye(:,2),'o');
title ('Restricted three body problem')
ylabel ('y(t)')
xlabel ('x(t)')
% --------------------------------------------------------------
function dydt = f(t,y)
mu = 1 / 82.45;
mustar = 1 - mu;
r13 = ((y(1) + mu)^2 + y(2)^2) ^ 1.5;
r23 = ((y(1) - mustar)^2 + y(2)^2) ^ 1.5;
dydt = [ y(3)
y(4)
2*y(4) + y(1) - mustar*((y(1)+mu)/r13) - ...
mu*((y(1)-mustar)/r23)
-2*y(3) + y(2) - mustar*(y(2)/r13) - mu*(y(2)/r23) ];
% --------------------------------------------------------------
function [value,isterminal,direction] = events(t,y)
% Locate the time when the object returns closest to the
% initial point y0 and starts to move away, and stop integration.
% Also locate the time when the object is farthest from the
% initial point y0 and starts to move closer.
%
% The current distance of the body is
%
% DSQ = (y(1)-y0(1))^2 + (y(2)-y0(2))^2
% = <y(1:2)-y0,y(1:2)-y0>
%
% A local minimum of DSQ occurs when d/dt DSQ crosses zero
% heading in the positive direction. We can compute d(DSQ)/dt as
%
% d(DSQ)/dt = 2*(y(1:2)-y0)'*dy(1:2)/dt = 2*(y(1:2)-y0)'*y(3:4)
%
y0 = [1.2; 0];
dDSQdt = 2 * ((y(1:2)-y0)' * y(3:4));
value = [dDSQdt; dDSQdt];
isterminal = [1; 0]; % Stop at local minimum
direction = [1; -1]; % [local minimum, local maximum]
例題：微分代数方程式
hb1daeは、hb1odeデモを微分代数方程式(DAE)問題として作り替えたものです。hb1ode
デモの中でコード化されている Robertson 問題は、スティッフな ODE
問題を解くコードに対する古典的なテスト問題です。
| 注意 Robertson 問題は、LSODI [3]へのプロローグへの例題です。 |
hb1odeで、問題は、初期条件
, 
, 
をもち、定常状態で解かれます。これらの微分方程式は、線形保存則を満足し。DAE
として問題を再定式化するのに使います。
明らかに、これらの方程式は、和が1にならない要素をもつ
に対する解をもっていません。
をもつ
の型に問題は、あります。
は明らかに正則でなく、hb1dae
は、ソルバに伝えません。ソルバは、問題が ODE ではなく、DAE
であることを識別します。同様に、整合性をもつ初期条件は明らかになりますが、例題では、
使い、整合性のある初期条件の計算を示しています。
MATLAB Help
ブラウザからこの例題を実行するには、デモ名をクリックするか、コマンドラインに、hb1dae と入力してください。hb1daeは、つぎの特徴をもっています。
により小さな絶対誤差のトレランスを課します。これは、
が他の要素より小さく、大きな変化が、相対的に短い時間で生じるためです。
の挙動をうまく示すために、付加的な点を加えます。
に
104を乗算して、
をプロットしたとき、見えるようにします。
function hb1dae
% HB1DAE スティッフな微分代数方程式(DAE)
% 定数、特異質量行列
M = [1 0 0
0 1 0
0 0 0];
% 初期化をテストするため、矛盾のある初期条件を使用
y0 = [1; 0; 1e-3];
tspan = [0 4*logspace(-6,6)];
% LSODI 例題の許容範囲を使用。'MassSingular' プロパティ
% は、DAE の自動検出をテストするデフォルト'maybe'のまま
% にします。
options = odeset('Mass',M,'RelTol',1e-4,...
'AbsTol',[1e-6 1e-10 1e-6],'Vectorized','on');
[t,y] = ode15s(@f,tspan,y0,options);
y(:,2) = 1e4*y(:,2);
semilogx(t,y);
ylabel('1e4 * y(:,2)');
title(['Robertson DAE problem with a Conservation Law, '...
'solved by ODE15S']);
xlabel('This is equivalent to the stiff ODEs coded in HB1ODE.');
% --------------------------------------------------------------
function out = f(t,y)
out = [ -0.04*y(1,:) + 1e4*y(2,:).*y(3,:)
0.04*y(1,:) - 1e4*y(2,:).*y(3,:) - 3e7*y(2,:).^2
y(1,:) + y(2,:) + y(3,:) - 1 ];
| |
ODE ソルバ性能の改良 | 質問と回答とトラブルシュート | |
