最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

ODE ソルバーの選択

常微分方程式

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

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

y''=9y

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

ODE のタイプ

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

  • y'=f(t,y) 形式の陽的 ODE。

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

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

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

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

関数 odeset を使用してオプション構造体を作成すると、一部のタイプの問題ではソルバーに追加情報を指定することもできます。

ODE 系

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

(y'1y'2y'n)=(f1(t,y1,y2,...,yn)f2(t,y1,y2,...,yn)fn(t,y1,y2,...,yn)),

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

{y'1=y2y'2=y1y22.

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

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

高次 ODE

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

y1=yy2=y'y3=y''yn=y(n1).

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

{y'1=y2y'2=y3y'n=f(t,y1,y2,...,yn).

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

y'''y''y+1=0.

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

y1=yy2=y'y3=y''

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

{y'1=y2y'2=y3y'3=y1y31.

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

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

複素数 ODE

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

y'=f(t,y),

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

yv=[Real(y)Imag(y)]fv=[Real(f(t,y))Imag(f(t,y))].

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

function f = complexf(t,y)
% Define function that takes and returns complex values
f = y.*t + 2*i;

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

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

解を取得するためのソルバーを実行すると、初期条件 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 および ode113ode45 よりも効率的なことがあります。

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

次の表に、さまざまなソルバーをそれぞれどのような場合に使用すればよいかについての一般的なガイダンスを示します。

ソルバー問題のタイプ精度使用時
ode45ノンスティッフ

ほとんどの場合に使用します。まず ode45 をソルバーとして使用してみてください。

ode23

粗い許容誤差が設定された問題や中程度のスティッフがある問題では、ode45 よりも ode23 の方が効率的なことがあります。

ode113低から高

厳しい許容誤差が設定された問題や ODE 関数の計算に時間がかかるような問題では、ode45 よりも ode113 の方が効率的なことがあります。

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 のサンプル ファイルと、それらで使用されているソルバーおよびオプションのリストを示します。例のサブセットがドキュメンテーションでも直接公開されている場合は、リンクを記載しています。

サンプル ファイル使用されているソルバーオプション指定説明ドキュメンテーション リンク
amp1daeode23t
  • 'Mass'

スティッフな DAE — 定数の特異質量行列を使用した電気回路

スティッフな微分代数方程式を解く
ballodeode23
  • 'Events'

  • 'OutputFcn'

  • 'OutputSel'

  • 'Refine'

  • 'InitialStep'

  • 'MaxStep'

簡単なイベント検出問題 — 跳ねるボール

ODE のイベント検出
batonodeode45
  • 'Mass'

時間と状態に依存した質量行列をもつ ODE — バトンの運動

brussodeode15s
  • 'JPattern'

  • 'Vectorized'

スティッフな大規模問題 — 化学反応の拡散作用 (Brusselator)

スティッフ ODE の求解
burgersodeode15s
  • 'Mass'

  • 'MStateDependence'

  • 'JPattern'

  • 'MvPattern'

  • 'RelTol'

  • 'AbsTol'

状態量に強く依存した質量行列をもつ ODE — 移動メッシュ法を使って Burger の方程式を解く

fem1odeode15s
  • 'Mass'

  • 'MStateDependence'

  • 'Jacobian'

時変の質量行列をもつスティッフな問題 — 有限要素法

fem2odeode23s
  • 'Mass'

定質量行列をもつスティッフな問題 — 有限要素法

hb1odeode15s

非常に長い区間で解くスティッフな ODE 問題 — Robertson 化学反応

hb1daeode15s
  • 'Mass'

  • 'RelTol'

  • 'AbsTol'

  • 'Vectorized'

保存則下でのスティッフで線形陰的な DAE — Robertson 化学反応

Robertson 問題を半陽的な微分代数方程式 (DAE) として求解
ihb1daeode15i
  • 'RelTol'

  • 'AbsTol'

  • 'Jacobian'

スティッフで完全陰的な DAE — Robertson 化学反応

ロバートソンの問題を陰的微分代数方程式 (DAE) として求解
iburgersodeode15i
  • 'RelTol'

  • 'AbsTol'

  • 'Jacobian'

  • 'JPattern'

陰的 ODE 系 — Burger の方程式

kneeodeode15s
  • 'NonNegative'

非負の制約をもつ "Knee 問題"

非負の ODE の解
orbitodeode45
  • 'RelTol'

  • 'AbsTol'

  • 'Events'

  • 'OutputFcn'

高度なイベント検出問題 — 制約付きの 3 体問題

ODE のイベント検出
rigidodeode45

ノンスティッフな問題 — 外力のない剛体のオイラー方程式

ノンスティッフ ODE の求解
vdpodeode15s
  • '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.

参考

|

関連するトピック

外部の Web サイト