ODE ソルバーの選択

常微分方程式

"常微分方程式" (ODE) には、1 つの独立変数 t に対する従属変数 y の導関数が 1 つ以上含まれます。t は通常、時間を表します。y の t に対す導関数を表すためにここで使用される表記法は、1 次導関数は $$y\text{'}$$、2 次導関数は $$y\text{'}\text{'}$$ などのようになります。ODE の "次数" は方程式における y の最高次の導関数に等しくなります。

たとえば、これは 2 次 ODE です。

$$y\text{'}\text{'}=9y$$

"初期値問題" で ODE の解を求めるには初期状態から開始します。初期条件 $$y_0$$ および解の取得対象期間 $$\left(t_0,t_f\right)$$ を使用して、解を反復的に求めます。各ステップで、ソルバーは前のステップの結果に対して特定のアルゴリズムを適用します。このような最初のステップでは、積分を進めるために必要な情報が初期状態によって提供されます。最終的な結果として、ODE ソルバーはタイム ステップのベクトル $$t=\left[t_0,t_1,t_2,\mathrm{...},t_f\right]$$ と、各ステップでの対応する解 $$y=\left[y_0,y_1,y_2,\mathrm{...},y_f\right]$$ を返します。

ODE のタイプ

MATLAB^® の ODE ソルバーでは次のタイプの 1 次 ODE の解を求めます。

• $$y\text{'}=f\left(t,y\right)$$ 形式の陽的 ODE。

• $$M\left(t,y\right)y\text{'}=f\left(t,y\right)$$ 形式の線形陰的な ODE。ここで、$$M\left(t,y\right)$$ は正則な質量行列です。この質量行列は時間または状態に依存するか、定数行列になります。線形陰的な ODE では、y の 1 次導関数の線形結合が使用されます。これらは質量行列にエンコードされます。

  線形陰的な ODE は陽的な形式 $$y\text{'}=M^{-1}\left(t,y\right)f\left(t,y\right)$$ に常に変換できます。ただし、ODE ソルバーに対して質量行列を直接指定すると、不便で計算量の多いこの変換を回避できます。

• $$y\text{'}$$ の一部の要素が欠けている場合、この方程式は "微分代数方程式" (DAE) と呼ばれます。DAE 系には "代数変数" がいくつか含まれます。代数変数は従属変数で、方程式にはその導関数がありません。DAE 系を等価な 1 次 ODE 系として書き換えるには、この方程式を微分して代数変数を除外します。DAE を ODE として書き換えるために必要な導関数の数は微分指数と呼ばれます。ode15s ソルバーと ode23t ソルバーでは指数が 1 の DAE の解を求めることができます。

• $$f\left(t,y,y\text{'}\right)=0$$ 形式の完全陰的な ODE。完全陰的な ODE は陽的な形式に書き換えることができません。また、代数変数がいくつか含まれることがあります。ode15i ソルバーは完全陰的な問題を対象に設計されています。指数が 1 の DAE もその対象に含まれます。

ODE 系

解を求める連立 ODE の方程式の数はいくつでも指定できます。原則として、方程式の数は利用できるコンピューター メモリのみによって制限を受けます。方程式系に n 個の方程式があるとします。

$$\left(\begin{array}{c}y{\text{'}}_1\\ y{\text{'}}_2\\ ⋮\\ y{\text{'}}_n\end{array}\right)=\left(\begin{array}{c}{f}_1\left(t,y_1,y_2,\mathrm{...},y_n\right)\\ f_2\left(t,y_1,y_2,\mathrm{...},y_n\right)\\ ⋮\\ f_n\left(t,y_1,y_2,\mathrm{...},y_n\right)\end{array}\right),$$

この場合、方程式をエンコードする関数は n 個の要素をもつベクトルを返し、それは $$y{\text{'}}_1,y{\text{'}}_2,\dots ,y{\text{'}}_n$$ の値に対応します。たとえば、次の 2 つを含む方程式系について考えます。

$$\{\begin{array}{l}y{\text{'}}_1=y_2\\ y{\text{'}}_2=y_1{y}_2-2.\end{array}$$

これらの方程式をエンコードする関数は次のようになります。

function dy = myODE(t,y)
  dy(1) = y(2);
  dy(2) = y(1)*y(2)-2;
end

高次 ODE

MATLAB ODE ソルバーでは 1 次方程式の解のみを求めることができます。高次 ODE は汎用の置き換え方法を使用して等価な 1 次の方程式系に書き換えなければなりません。

$$\begin{array}{l}{y}_1=y\\ y_2=y\text{'}\\ y_3=y\text{'}\text{'}\\ ⋮\\ y_n=y^{(n-1)}.\end{array}$$

これらの置き換えの結果は、次のような n 個の 1 次方程式系になります。

$$\{\begin{array}{l}y{\text{'}}_1=y_2\\ y{\text{'}}_2=y_3\\ ⋮\\ y{\text{'}}_n=f\left(t,y_1,y_2,\mathrm{...},y_n\right).\end{array}$$

たとえば、次のような 3 次 ODE について考えます。

$$y\text{'}\text{'}\text{'}-y\text{'}\text{'}y+1=0.$$

次のような置き換えを使用します。

$$\begin{array}{l}{y}_1=y\\ y_2=y\text{'}\\ y_3=y\text{'}\text{'}\end{array}$$

この結果、等価な 1 次方程式系は次のようになります。

$$\{\begin{array}{l}y{\text{'}}_1=y_2\\ y{\text{'}}_2=y_3\\ y{\text{'}}_3=y_1{y}_3-1.\end{array}$$

この方程式系のコードは次のようになります。

function dydt = f(t,y)
  dydt(1) = y(2);
  dydt(2) = y(3);
  dydt(3) = y(1)*y(3)-1;
end

複素数 ODE

次の複素数 ODE 方程式について考えます。

$$y\text{'}=f\left(t,y\right),$$

ここで $$y=y_1+i{y}_2$$ です。この解を求めるには、実数部と虚数部を別々の解の要素に分けてから、最後にその結果を再結合します。概念的には次のようになります。

$$\begin{array}{l}yv=\left[\text{Real}\left(y\right)\text{Imag}\left(y\right)\right]\\ fv=\left[\text{Real}\left(f\left(t,y\right)\right)\text{Imag}\left(f\left(t,y\right)\right)\right].\end{array}$$

たとえば、ODE が $$y\text{'}=yt+2i$$ の場合、関数ファイルを使用してこの方程式を表すことができます。

function f = complexf(t,y)
  f = y.*t + 2*i;
end

この場合、実数部と虚数部を分けるコードは次のようになります。

function fv = imaginaryODE(t,yv)
  % Construct y from the real and imaginary components
  y = yv(1) + i*yv(2);            

  % Evaluate the function
  yp = complexf(t,y);             

  % Return real and imaginary in separate components
  fv = [real(yp); imag(yp)]; 
end     

解を取得するためのソルバーを実行すると、初期条件 y0 も実数部と虚数部に分けられ、解の要素それぞれに対して初期条件が提供されます。

y0 = 1+i;
yv0 = [real(y0); imag(y0)];
tspan = [0 2];
[t,yv] = ode45(@imaginaryODE, tspan, yv0);

解を取得したら、実数部と虚数部の両方の要素を結合して最終結果を得ます。

y = yv(:,1) + i*yv(:,2);

基本的なソルバー選択

ODE のほとんどの問題には ode45 の実行が適しているため、通常はまずこのソルバーを選択してください。ただし、精度に関する要件が緩い場合または厳しい場合は、ode23、ode78、ode89、および ode113 が ode45 よりも効率的なことがあります。

ODE の一部の問題では、"スティッフ性" が見られる場合、つまり、評価が困難な場合があります。スティッフとは明確な定義が難しい用語ですが、一般的には、問題中のいずれかの部分でスケールに差があるときにスティッフになります。たとえば、時間スケールに大幅な違いがあると変化する解の要素が ODE に 2 つある場合、その方程式はスティッフな可能性があります。ノンスティッフ ソルバー (ode45 など) で問題の解を求められない場合や極端に処理が遅い場合は、その問題はスティッフであると判定できます。ノンスティッフ ソルバーの処理が非常に遅いことがわかった場合は、代わりに ode15s などのスティッフ ソルバーを使用してみてください。スティッフ ソルバーを使用する場合、ヤコビ行列やそのスパース パターンを渡すことによって信頼性と効率性を高めることができます。

ode オブジェクトを使用して、問題のプロパティに基づいてソルバーの選択を自動化できます。どのソルバーを使用すればよいかわからない場合のために、下表に各ソルバーをどのような場合に使用すればよいかについての一般的なガイドラインを示します。

ソルバー	問題のタイプ	精度	使用時
ode45	ノンスティッフ	中	ほとんどの場合に使用します。まず ode45 をソルバーとして使用してみてください。
ode23		低	粗い許容誤差が設定された問題や中程度のスティッフがある問題では、ode45 よりも ode23 の方が効率的なことがあります。
ode113		低から高	厳しい許容誤差が設定された問題や ODE 関数の計算に時間がかかるような問題では、ode45 よりも ode113 の方が効率的なことがあります。
ode78		高	精度に関する要件が高い滑らかな解の問題では、ode45 よりも ode78 の方が効率的なことがあります。
ode89		高	積分の時間間隔が長い場合や許容誤差が特に厳しい場合など、非常に滑らかな問題では、ode78 よりも ode89 の方が効率的なことがあります。
ode15s	スティッフ	低から中	ode45 が失敗する場合や非効率な場合にその問題がスティッフであることが疑われるときは、ode15s を試してください。ode15s は微分代数方程式 (DAE) の解を求める場合にも使用します。
ode23s		低	粗い許容誤差が設定された問題では、ode15s よりも ode23s の方が効率的なことがあります。これによって、ode15s が効率的ではないときに、スティッフな問題の解を求めることができる場合があります。 ode23s では各ステップでヤコビアンを計算します。そのため、odeset を介してヤコビアンを渡すと、効率性と精度を最大化できるため、高い効果が得られます。 質量行列がある場合、これは定数でなければなりません。
ode23t		低	ode23t は、問題のスティッフが中程度に留まり、数値的減衰のない解が必要な場合に使用します。 ode23t では微分代数方程式 (DAE) の解を求めることができます。
ode23tb		低	ode23s と同様に、粗い許容誤差が設定された問題では、ode15s よりも ode23tb ソルバーの方が効率的なことがあります。
ode15i	完全陰的	低	完全陰的な問題 f(t,y,y’) = 0 の場合や指数 1 の微分代数方程式 (DAE) の場合は、ode15i を使用します。

各ソルバーをどのような場合に使用すればよいかの詳細と追加推奨事項については、[5]を参照してください。

ODE の例とファイルの概要

ODE のほとんどの問題に対して開始点として効果的に使用できるサンプル ファイルがいくつかあります。[微分方程式の例] アプリでは例を簡単に参照したり実行したりできます。このアプリを実行するには、次のように入力します。

odeexamples

個々のサンプル ファイルを編集用に開くには、次のように入力します。

edit exampleFileName.m

例を実行するには、次のように入力します。

exampleFileName

次の表に、利用可能な ODE および DAE のサンプル ファイルと、それらで使用されているソルバーおよびオプションのリストを示します。例のサブセットがドキュメンテーションでも直接公開されている場合は、リンクを記載しています。

サンプル ファイル	使用されているソルバー	オプション指定	説明	ドキュメンテーション リンク
amp1dae	ode23t	'Mass'	スティッフな DAE — 定数の特異質量行列を使用した電気回路	トランジスタのスティッフな微分代数方程式の求解
ballode	ode23	'Events' 'OutputFcn' 'OutputSel' 'Refine' 'InitialStep' 'MaxStep'	簡単なイベント検出問題 — 跳ねるボール	ODE のイベント検出
batonode	ode45	'Mass'	時間と状態に依存した質量行列をもつ ODE — バトンの運動	空中に投げられたバトンの動きを表す方程式の求解
brussode	ode15s	'JPattern' 'Vectorized'	スティッフな大規模問題 — 化学反応の拡散作用 (Brusselator)	スティッフ ODE の求解
burgersode	ode15s	'Mass' 'MStateDependence' 'JPattern' 'MvPattern' 'RelTol' 'AbsTol'	状態量に強く依存した質量行列をもつ ODE — 移動メッシュ法を使って Burger の方程式を解く	状態依存の強い質量行列による ODE の求解
fem1ode	ode15s	'Mass' 'MStateDependence' 'Jacobian'	時変の質量行列をもつスティッフな問題 — 有限要素法	—
fem2ode	ode23s	'Mass'	定質量行列をもつスティッフな問題 — 有限要素法	—
hb1ode	ode15s	—	非常に長い区間で解くスティッフな ODE 問題 — Robertson 化学反応	—
hb1dae	ode15s	'Mass' 'RelTol' 'AbsTol' 'Vectorized'	保存則下でのスティッフで線形陰的な DAE — Robertson 化学反応	Robertson 問題を半陽的な微分代数方程式 (DAE) として求解
ihb1dae	ode15i	'RelTol' 'AbsTol' 'Jacobian'	スティッフで完全陰的な DAE — Robertson 化学反応	ロバートソンの問題を陰的微分代数方程式 (DAE) として求解
iburgersode	ode15i	'RelTol' 'AbsTol' 'Jacobian' 'JPattern'	陰的 ODE 系 — Burger の方程式	—
kneeode	ode15s	'NonNegative'	非負の制約をもつ "Knee 問題"	非負の ODE の解
orbitode	ode45	'RelTol' 'AbsTol' 'Events' 'OutputFcn'	高度なイベント検出問題 — 制約付きの 3 体問題	ODE のイベント検出
rigidode	ode45	—	ノンスティッフな問題 — 外力のない剛体のオイラー方程式	ノンスティッフ ODE の求解
vdpode	ode15s	'Jacobian'	パラメーター化したファン デル ポールの方程式 (μ が大きいときスティッフ)	スティッフ ODE の求解

参照

[1] Shampine, L. F. and M. K. Gordon, Computer Solution of Ordinary Differential Equations: the Initial Value Problem, W. H. Freeman, San Francisco, 1975.

[2] Forsythe, G., M. Malcolm, and C. Moler, Computer Methods for Mathematical Computations, Prentice-Hall, New Jersey, 1977.

[3] Kahaner, D., C. Moler, and S. Nash, Numerical Methods and Software, Prentice-Hall, New Jersey, 1989.

[4] Shampine, L. F., Numerical Solution of Ordinary Differential Equations, Chapman & Hall, New York, 1994.

[5] Shampine, L. F. and M. W. Reichelt, “The MATLAB ODE Suite,” SIAM Journal on Scientific Computing, Vol. 18, 1997, pp. 1–22.

[6] Shampine, L. F., Gladwell, I. and S. Thompson, Solving ODEs with MATLAB, Cambridge University Press, Cambridge UK, 2003.

参考

ode | odeset | odextend

関連するトピック

• ノンスティッフ ODE の求解
• スティッフ ODE の求解
• 微分代数方程式 (DAE) の求解

外部の Web サイト

• 常微分方程式