ode15i

完全陰的微分方程式の求解 — 可変次数法

構文

[t,y] = ode15i(odefun,tspan,y0,yp0)

[t,y] = ode15i(odefun,tspan,y0,yp0,options)

[t,y,te,ye,ie] = ode15i(odefun,tspan,y0,yp0,options)

sol = ode15i(___)

説明

[t,y] = ode15i(odefun,tspan,y0,yp0) は、tspan = [t0 tf] のときに、初期条件 y0 および yp0 を使用して、連立微分方程式 $$f\left(t,y,y\text{'}\right)=0$$ を t0 から tf まで積分します。解の配列 y の各行は、列ベクトル t に返される値に対応します。

[t,y] = ode15i(odefun,tspan,y0,yp0,options) は options (関数 odeset を使用して作成された引数) で定義された積分設定も使用します。たとえば、AbsTol オプションおよび RelTol オプションを使用して絶対許容誤差と相対許容誤差を指定したり、Jacobian オプションを使用してヤコビ行列を指定することができます。

[t,y,te,ye,ie] = ode15i(odefun,tspan,y0,yp0,options) はさらに、(t,y,y') の関数 (イベント関数) がゼロになる点を求めます。出力の te はイベント時点、ye はイベント時点における解、ie はトリガーされたイベントのインデックスです。

各関数に対して、ゼロで積分を終了するかどうかと、ゼロクロッシングの方向を考慮するかどうかを指定します。これを行うには、myEventFcn や @myEventFcn などの関数に 'Events' プロパティを設定し、対応する関数 [value,isterminal,direction] = myEventFcn(t,y,yp) を作成します。詳細については、ODE のイベント検出を参照してください。

sol = ode15i(___) は、区間 [t0 tf] の任意の点で解を計算する関数 deval で使用できる構造体を返します。前述の構文にある任意の入力引数の組み合わせが使用できます。

例

Weissinger の陰的 ODE の求解

矛盾のない初期条件を計算し、陰的 ODE を ode15i で解きます。

Weissinger の方程式は次のとおりです。

${\mathit{t}\mathit{y}}^2{\left({\mathit{y}}^{\prime }\right)}^3-{\mathit{y}}^3{\left({\mathit{y}}^{\prime }\right)}^2+\mathit{t}\left({\mathit{t}}^2+1\right){\mathit{y}}^{\prime }-{\mathit{t}}^2\mathit{y}=0$.

方程式が一般的な形式 $\mathit{f}\left(\mathit{t},\mathit{y},{\mathit{y}}^{\prime }\right)=0$ であるため、関数 ode15i を使用して陰的微分方程式を解くことができます。

方程式のコード化

ode15i に適した形式で方程式のコードを作成するには、$\mathit{t}$、$\mathit{y}$、および ${\mathit{y}}^{\prime }$ への入力を含み、方程式の残差値を返す関数を記述する必要があります。関数 @weissinger でこの方程式をエンコードします。関数ファイルを表示します。

type weissinger

function res = weissinger(t,y,yp)
%WEISSINGER  Evaluate the residual of the Weissinger implicit ODE
%
%   See also ODE15I.

%   Jacek Kierzenka and Lawrence F. Shampine
%   Copyright 1984-2014 The MathWorks, Inc.

res = t*y^2 * yp^3 - y^3 * yp^2 + t*(t^2 + 1)*yp - t^2 * y;

矛盾のない初期条件の計算

ode15i ソルバーは "矛盾のない" 初期条件を必要とします。つまり、ソルバーに提供される初期条件は以下を満たさなければなりません。

$\mathit{f}\left({\mathit{t}}_0,\mathit{y},{\mathit{y}}^{\prime }\right)=0$.

矛盾した初期条件を提供することが可能であり、また ode15i は整合性を確認しないため、補助関数 decic を使用してこのような条件を計算することを推奨します。decic は指定されたいくつかの変数を固定し、固定されていない変数の矛盾のない初期条件を計算します。

ここでは、初期値 $\mathit{y}\left({\mathit{t}}_0\right)=\sqrt{\frac{3}{2}}$ を固定し、${\mathit{y}}^{\prime }\left({\mathit{t}}_0\right)=0$ の初期推定から計算を開始して、decic で導関数 ${\mathit{y}}^{\prime }\left({\mathit{t}}_0\right)$ の矛盾のない初期条件を計算します。

t0 = 1;
y0 = sqrt(3/2);
yp0 = 0;
[y0,yp0] = decic(@weissinger,t0,y0,1,yp0,0)

y0 = 
1.2247

yp0 = 
0.8165

方程式の求解

decic によって返された矛盾のない初期条件を ode15i で使用して、ODE を時間区間 $\left[110\right]$ について解きます。

[t,y] = ode15i(@weissinger,[1 10],y0,yp0);

結果のプロット

この ODE の厳密解は次のとおりです。

$\mathit{y}\left(\mathit{t}\right)=\sqrt{{\mathit{t}}^2+\frac{1}{2}}$.

ode15i によって計算される数値解 y を解析解 ytrue に対してプロットします。

ytrue = sqrt(t.^2 + 0.5);
plot(t,y,'*',t,ytrue,'-o')
legend('ode15i', 'exact')

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

この例では、ある ODE 系を完全陰的微分代数方程式 (DAE) として再定式化します。hb1ode.m でコード化されているロバートソンの問題は、スティッフな ODE 問題を解くプログラム向けの古典的なテスト問題です。連立方程式は次のとおりです。

hb1ode では、初期条件 、、および を使用して、この ODE 系の定常状態の解を求めます。ただし、これらの方程式は次の線形保存則も満たします。

この解と初期条件に基づくと、保存則は次のようになります。

この保存則を使用して の状態を決定することにより、この問題を DAE 系として書き換えることができます。これにより、問題が次の陰的 DAE 系として再定式化されます。

関数 robertsidae でこの DAE 系をエンコードします。

function res = robertsidae(t,y,yp)
res = [yp(1) + 0.04*y(1) - 1e4*y(2)*y(3);
   yp(2) - 0.04*y(1) + 1e4*y(2)*y(3) + 3e7*y(2)^2;
   y(1) + y(2) + y(3) - 1];

ロバートソンの問題のこの定式化を含む完全なサンプル コードは ihb1dae.m として提供されています。

許容誤差と の値を設定します。

options = odeset('RelTol',1e-4,'AbsTol',[1e-6 1e-10 1e-6], ...
   'Jacobian',{[],[1 0 0; 0 1 0; 0 0 0]});

decic を使用して推定値から矛盾のない初期条件を計算します。y0 の最初の 2 つの要素を固定して、この問題を半陽的 DAE 系として定式化した hb1dae.m で ode15s が求めたものと同一の矛盾のない初期条件を求めます。

y0 = [1; 0; 1e-3];
yp0 = [0; 0; 0];
[y0,yp0] = decic(@robertsidae,0,y0,[1 1 0],yp0,[],options);

ode15i を使用して DAE 系を求解します。

tspan = [0 4*logspace(-6,6)];
[t,y] = ode15i(@robertsidae,tspan,y0,yp0,options);

解の要素をプロットします。2 番目の解の要素が他と比べて小さいため、プロットする前に 1e4 を乗算します。

y(:,2) = 1e4*y(:,2);
semilogx(t,y)
ylabel('1e4 * y(:,2)')
title('Robertson DAE problem with a Conservation Law, solved by ODE15I')

入力引数

odefun — 解を求める関数

関数ハンドル

解を求める関数。積分する関数を定義する関数ハンドルとして指定します。

スカラー t、列ベクトル y および yp をとる関数 f = odefun(t,y,yp) は、データ型が single または double で $$f\left(t,y,y\text{'}\right)$$ に対応する列ベクトル f を返さなければなりません。odefun は、t、y、yp のいずれかの入力が関数に使用されない場合でも、それら 3 つの入力を受け入れなければなりません。

たとえば、$$y\text{'}-y=0$$ を解くには、次の関数を使用します。

function f = odefun(t,y,yp)
  f = yp - y;
end

連立方程式の場合、odefun の出力はベクトルです。ベクトルの各要素は 1 つの方程式の計算値です。たとえば、次の 2 つを含む連立方程式について考えます。

$$\begin{array}{l}y{\text{'}}_1-y_2=0\\ y{\text{'}}_2+1=0,\end{array}$$

各タイム ステップで各方程式の値を計算する関数は以下のとおりです。

function dy = odefun(t,y,yp)
  dy = zeros(2,1);
  dy(1) = yp(1)-y(2);
  dy(2) = yp(2)+1;
end

関数 odefun に追加パラメーターを指定する方法の詳細については、関数のパラメーター化を参照してください。

例: @myFcn

データ型: function_handle

tspan — 積分区間

ベクトル

積分区間。ベクトルとして指定します。少なくとも、tspan は開始時点と終了時点を指定する 2 要素ベクトル [t0 tf] でなければなりません。t0 と tf の間の特定時点における解を取得するには、[t0,t1,t2,...,tf] の形式の長いベクトルを使用します。tspan の要素は、単調増加または単調減少でなければなりません。

ソルバーは開始時点 tspan(1) で y0 によって指定される初期条件を設定し、tspan(1) から tspan(end) まで積分します。

• tspan に 2 つの要素がある場合 ([t0 tf])、ソルバーは区間内の個々の内部積分ステップで計算した解を返します。

• tspan に 2 つを超える要素がある場合 ([t0,t1,t2,...,tf])、ソルバーは指定された各点で計算した解を返します。ただし、ソルバーは tspan で指定された各点に正確にステップするわけではありません。代わりに、ソルバーは独自の内部ステップを使用して解を計算し、tspan 内の要求された点で解を評価します。指定された点で出力された解の精度の次数は、各内部ステップで計算された解と同じです。

  中間点をいくつか指定しても計算効率にほとんど影響しませんが、大規模な方程式の場合はメモリ管理に影響を及ぼす可能性があります。

ソルバーは tspan の値を使用して InitialStep と MaxStep に適した値を計算します。

• tspan に複数の中間点が含まれる場合 ([t0,t1,t2,...,tf])、指定した点は問題のスケールの目安となり、これはソルバーが使用する InitialStep の値に影響することがあります。そのため、ソルバーで得られる解は、tspan を 2 要素ベクトルとして指定するか、中間点を含むベクトルとして指定するかによって異なる場合があります。

• tspan の最初と最後の値は、最大ステップ サイズ MaxStep の計算に使用されます。そのため、tspan の最初または最後の値を変更すると、ソルバーが別のステップ シーケンスを使用する可能性があり、そのため解が変わることがあります。

例: [1 10]

例: [1 3 5 7 9 10]

データ型: single | double

y0 — y の初期条件

ベクトル

y の初期条件。ベクトルとして指定します。odefun に定義された各方程式の初期条件を y0 に含めるために、y0 は odefun のベクトル出力と同じ長さでなければなりません。

y0 と yp0 の初期条件に矛盾があってはなりません。つまり $$f\left(t_0,y_0,y{\text{'}}_0\right)=0$$ です。矛盾がなく、推定値に近い初期条件を計算するには、関数 decic を使用します。

データ型: single | double

yp0 — y’ の初期条件

列ベクトル

y’ の初期条件。列ベクトルとして指定します。odefun に定義された各変数の初期条件を yp0 に含めるために、yp0 は odefun のベクトル出力と同じ長さでなければなりません。

y0 と yp0 の初期条件に矛盾があってはなりません。つまり $$f\left(t_0,y_0,y{\text{'}}_0\right)=0$$ です。矛盾がなく、推定値に近い初期条件を計算するには、関数 decic を使用します。

データ型: single | double

options — オプション構造体

構造体配列

オプション構造体。構造体配列として指定します。オプション構造体の作成または変更を行うには、関数 odeset を使用します。

各 ODE ソルバーと互換性のあるオプションの一覧については、ODE オプションの概要を参照してください。

例: options = odeset('RelTol',1e-5,'Stats','on','OutputFcn',@odeplot) は、相対許容誤差 1e-5 を指定し、ソルバー統計の表示をオンにして、解の計算時に解をプロットする出力関数 @odeplot を指定します。

データ型: struct

出力引数

t — 評価点

列ベクトル

評価点。列ベクトルとして返されます。

• tspan に 2 つの要素がある場合 ([t0 tf])、t には積分の実行に使用される内部評価点が含まれます。

• tspan に 2 つより多い要素が含まれる場合、t は tspan と同じです。

y — 解

配列

解。配列として返されます。y の各行は、t の対応する行に返される値での解です。

te — イベント時点

列ベクトル

イベント時点。列ベクトルとして返されます。te のイベント時点は ye に返された解に対応し、ie は発生したイベントを指定します。

ye — イベント時点での解

配列

イベント時点での解。配列として返されます。te のイベント時点は ye に返された解に対応し、ie は発生したイベントを指定します。

ie — トリガーされるイベント関数のインデックス

列ベクトル

トリガーされるイベント関数のインデックス。列ベクトルとして返されます。te のイベント時点は ye に返された解に対応し、ie は発生したイベントを指定します。

sol — 評価対象の構造体

構造体配列

評価対象の構造体。構造体配列として返されます。この構造体を関数 deval と共に使用して区間 [t0 tf] の任意の点で解を評価します。構造体配列 sol は、常に以下のフィールドを含みます。

構造体フィールド	説明
sol.x	ソルバーによって選択されたステップの行ベクトル。
sol.y	解。各列の sol.y(:,i) には、時点 sol.x(i) での解が含まれます。
sol.solver	ソルバー名。

さらに、odeset の Events オプションを指定してイベントが検出された場合、sol は以下のフィールドも含みます。

構造体フィールド	説明
sol.xe	イベントが発生した点。sol.xe(end) には、終了イベントの厳密な点 (存在する場合) が含まれます。
sol.ye	sol.xe のイベントに対応する解。
sol.ie	Events オプションに指定された関数により返されるベクトルのインデックス。値は、どのイベントをソルバーが検出したかを示します。

ヒント

• ode15i にヤコビ行列を指定することは、信頼性と効率性のために重要です。また、方程式が大規模かつスパースである場合、ヤコビ スパース パターンを指定してソルバーを補助することもできます。いずれの場合でも、odeset と Jacobian オプションまたは JPattern オプションを使用して行列を渡します。

アルゴリズム

ode15i は 1 次～ 5 次の後退差分式 (BDF) に基づく可変ステップ、可変次数 (VSVO) ソルバーです。ode15i は完全陰的微分方程式およびインデックス 1 の微分代数方程式 (DAE) と使用できるように設計されています。補助関数 decic は ode15i [1] との使用に適した矛盾のない初期条件を計算します。

参照

[1] Lawrence F. Shampine, “Solving 0 = F(t, y(t), y′(t)) in MATLAB,” Journal of Numerical Mathematics, Vol.10, No.4, 2002, pp. 291-310.

バージョン履歴

R2006a より前に導入

参考

decic | ode15s | ode23t | odeset | odeget | deval

トピック

• ODE ソルバーの選択
• ODE オプションの概要
• 微分代数方程式 (DAE) の求解
• 無名関数
• ODE の一般的な問題のトラブルシューティング