ドキュメンテーション

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

decic

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

構文

  • [y0_new,yp0_new] = decic(odefun,t0,y0,fixed_y0,yp0,fixed_yp0)
  • [y0_new,yp0_new] = decic(odefun,t0,y0,fixed_y0,yp0,fixed_yp0,options)
  • [y0_new,yp0_new,resnrm] = decic(___)

説明

[y0_new,yp0_new] = decic(odefun,t0,y0,fixed_y0,yp0,fixed_yp0) は完全陰関数 odefun の初期条件の推定値として y0yp0 を使用し、fixed_y0 および fixed_yp0 で指定された成分を固定して、固定されていない成分の値を計算します。結果は、矛盾のない初期条件の完全なセットです。新しい値 yo_new および yp0_newodefun(t0,y0_new,yp0_new) = 0 を満たし、ode15i の初期条件としての使用に適します。

[y0_new,yp0_new] = decic(odefun,t0,y0,fixed_y0,yp0,fixed_yp0,options) は、AbsTol および RelTol の値を指定する options 構造体 options も使用します。options 構造体を作成するには odeset を使用します。

[y0_new,yp0_new,resnrm] = decic(___) は、odefun(t0,y0_new,yp0_new) のノルムを resnrm として返します。ノルムが大きすぎると思われる場合は、options を使用して相対許容誤差 RelTol を小さくします (既定値は 1e-3 です)。

すべて折りたたむ

陰的な方程式系を考えます。

$$\begin{array}{cc} 0 &= 2y'_1 - y_2\\ 0 &=y_1+y_2 \end{array}$$

これらの方程式は十分わかりやすいため、変数の矛盾のない初期条件を簡単に読み取ることができます。たとえば、 $y_1 = 1$ を固定すると、2 番目の方程式から $y_2=-1$ であり、1 番目の方程式から $y'_1 = -1/2$ です。 $y_1$ $y'_1$ および $y_2$ のこれらの値は方程式を満たすので、矛盾しません。

decic を使用し、値 $y_1 = 1$ を固定して矛盾のない初期条件を計算し、これらの値を確認します。y0 = [1 0] および yp0 = [0 0] の推定を使用すると方程式は満たされず、したがって矛盾しています。

odefun = @(t,y,yp) [2*yp(1)-y(2); y(1)+y(2)];
t0 = 0;
y0 = [1 0];
yp0 = [0 0];
[y0,yp0] = decic(odefun,t0,y0,[1 0],yp0,[])
y0 =

     1
    -1


yp0 =

   -0.5000
         0

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

入力引数

すべて折りたたむ

求解する関数。積分する関数を定義する関数ハンドルとして指定します。odefun は、ode15i を使用して求解する陰的微分方程式系を表します。

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

たとえば、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

初期時刻。スカラーとして指定します。decic は初期時刻を使用して、odefun(t0,y0_new,yp0_new) = 0 を満たす矛盾のない初期条件を計算します。

データ型: single | double

y 成分の初期推定値。ベクトルとして指定します。y0 の各成分は、odefun で定義される方程式系の 1 つの従属変数 yn の初期条件を指定します。

データ型: single | double

固定する y 成分。1 と 0 からなるベクトル、あるいは [] として指定します。

  • y0(i) の推定値の変更が許可されない場合は fixed_y0(i) = 1 を設定します。

  • 任意のエントリが変更可能な場合は fixed_y0 = [] を設定します。

length(yp0) 個を超える成分を固定することはできません。問題によっては、y0 または yp0 の特定の成分を常に固定できるとは限りません。必要な数より多くの成分を固定しないことをお勧めします。

y' 成分の初期推定値。ベクトルとして指定します。yp0 の各成分は、odefun で定義される方程式系の 1 つの微分された従属変数 y'n の初期条件を指定します。

データ型: single | double

固定する y' 成分。1 と 0 からなるベクトル、あるいは [] として指定します。

  • yp0(i) の推定値の変更が許可されない場合は fixed_yp0(i) = 1 を設定します。

  • 任意のエントリが変更可能な場合は fixed_yp0 = [] を設定します。

length(yp0) 個を超える成分を固定することはできません。問題によっては、y0 または yp0 の特定の成分を常に固定できるとは限りません。必要な数より多くの成分を固定しないことをお勧めします。

オプション構造体。構造体配列として指定します。オプション構造体の作成または変更を行うには、関数 odeset を使用します。関数 decic で使用する関連オプションには RelTolAbsTol があり、これらは初期条件の計算に使用する誤差のしきい値を制御します。

例: options = odeset('RelTol',1e-5)

データ型: struct

出力引数

すべて折りたたむ

y0 の矛盾のない初期条件。ベクトルとして返されます。値 resnrm が小さい場合、yo_new および yp0_newodefun(t0,y0_new,yp0_new) = 0 を満たし、ode15i の初期条件としての使用に適します。

yp0 の矛盾のない初期条件。ベクトルとして返されます。値 resnrm が小さい場合、yo_new および yp0_newodefun(t0,y0_new,yp0_new) = 0 を満たし、ode15i の初期条件としての使用に適します。

残差のノルム。ベクトルとして返されます。resnrmodefun(t0,y0_new,yp0_new) のノルムです。

  • resnrm が小さい場合、odefun(t0,y0_new,yp0_new) = 0 を満たす矛盾のない初期条件が decic で正しく計算されたことを示します。

  • resnrm が大きい場合は、options の入力を使用して、誤差のしきい値 RelTol および AbsTol を調整してください。

詳細

すべて折りたたむ

ヒント

  • ihb1dae および iburgersode の例のファイルは、decic を使用して初期条件を計算してから ode15i を使用して求解します。コードを表示するには、edit ihb1dae または edit iburgersode を入力します。

  • さらに、decic を使用して、ode15s または ode23t により解を求める DAE の矛盾のない初期条件を計算することもできます。これを行うには、次の手順に従ってください。

    1. 方程式系を完全に陰的な形式 f(t,y,y') = 0 に書き換えます。

    2. decic を呼び出して、これらの方程式の矛盾のない初期条件を計算します。

    3. ソルバーの呼び出し内に初期条件として y0_new を指定し、odesetInitialSlope オプションの値として yp_new を指定します。

参考

| |

R2006a より前に導入

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