ドキュメンテーション

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

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(t,y,y')=0t0 から tf まで積分します。解の配列 y の各行は、列ベクトル t に返される値に対応します。

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

[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 で使用できる構造体を返します。前述の構文にある任意の入力引数の組み合わせが使用できます。

すべて折りたたむ

decic を使用して、Weissinger の陰的 ODE の矛盾のない初期条件を計算します。decicy(t0) の初期値を固定して、矛盾のない y'(t0) の初期値を計算します。関数 weissinger は、陰的 ODE の残差を評価します。

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

decic によって返された結果を ode15i で使用して ODE を解きます。数値解 y を解析解 ytrue に対してプロットします。

[t,y] = ode15i(@weissinger,[1 10],y0,yp0);
ytrue = sqrt(t.^2 + 0.5);
plot(t,y,t,ytrue,'o')

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

$$\begin{array}{cl} y'_1 &= -0.04y_1 + 10^4 y_2y_3\\ y'_2 &= 0.04y_1 -
10^4 4y_2y_3-(3 \times 10^7)y_2^2\\ y'_3 &= (3 \times
10^7)y_2^2.\end{array}$$

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

$$y'_1 + y'_2 + y'_3 = 0.$$

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

$$y_1 + y_2 + y_3 = 1.$$

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

$$\begin{array}{cl} 0 &= y'_1 +0.04y_1 - 10^4 y_2y_3\\ 0 &= y'_2 -0.04y_1
+ 10^4 y_2y_3+(3 \times 10^7)y_2^2\\ 0 &= y_1 + y_2 + y_3 -
1.\end{array}$$

関数 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 として提供されています。

許容誤差と $\partial f / \partial y'$ の値を設定します。

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.mode15s が求めたものと同一の矛盾のない初期条件を求めます。

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

入力引数

すべて折りたたむ

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

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

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

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

方程式系の場合、odefun の出力はベクトルです。各方程式が解ベクトルの 1 つの要素になります。たとえば、

y'1y2=0y'2+1=0,

を解くには、次の関数を使用します。

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

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

例: @myFcn

データ型: function_handle

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

ソルバーは、tspan(1) での初期条件 y0 を設定し、tspan(1) から tspan(end) まで積分します。

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

  • tspan に 2 個より多い要素 [t0,t1,t2,...,tf] が含まれる場合、ソルバーは指定された各点で計算した解を返します。これはソルバーが tspan(1) から tspan(end) に進むときに使用する内部ステップに影響 "しません"。したがって、ソルバーは必ずしも tspan で指定された各点に正確にステップするとは限りません。しかし、指定された点で出力された解の精度の次数は、各内部ステップで計算された解と同じです。

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

ソルバーで得られる解は、tspan を 2 要素ベクトルとして指定するか、中間点を含むベクトルとして指定するかによって異なる場合があります。tspan にいくつかの中間点が含まれる場合、これは問題のスケールを示し、ソルバーが使用する初期ステップのサイズに影響を及ぼすことがあります。

例: [1 10]

例: [1 3 5 7 9 10]

データ型: single | double

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

y0yp0 の初期条件に矛盾があってはなりません。つまり f(t0,y0,y'0)=0 です。矛盾がなく、推定値に近い初期条件を計算するには、関数 decic を使用します。

データ型: single | double

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

y0yp0 の初期条件に矛盾があってはなりません。つまり f(t0,y0,y'0)=0 です。矛盾がなく、推定値に近い初期条件を計算するには、関数 decic を使用します。

データ型: single | double

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

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

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

データ型: struct

出力引数

すべて折りたたむ

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

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

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

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

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

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

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

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

構造体フィールド説明

sol.x

ソルバーによって選択されたステップの行ベクトル。

sol.y

解。各列の sol.y(:,i) には、時点 sol.x(i) での解が含まれます。

sol.solver

ソルバー名。

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

構造体フィールド説明

sol.xe

イベントが発生した点。sol.xe(end) には、終了イベントの厳密な点 (存在する場合) が含まれます。

sol.ye

sol.xe のイベントに対応する解。

sol.ie

Events オプションに指定された関数により返されるベクトルのインデックス。値は、どのイベントをソルバーが検出したかを示します。

詳細

すべて折りたたむ

ヒント

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

アルゴリズム

ode15i は 1 次~ 5 次の後退差分式 (BDF) に基づく可変ステップ、可変次数 (VSVO) ソルバーです。ode15i は完全陰的微分方程式およびインデックス 1 の微分代数方程式 (DAE) と使用できるように設計されています。補助関数 decicode15i [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 より前に導入

この情報は役に立ちましたか?