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 の初期条件の推定値として y0 と yp0 を使用し、fixed_y0 および fixed_yp0 で指定された成分を固定して、固定されていない成分の値を計算します。結果は、矛盾のない初期条件の完全なセットです。新しい値 yo_new および yp0_new は odefun(t0,y0_new,yp0_new) = 0 を満たし、ode15i の初期条件としての使用に適します。

[y0_new,yp0_new] = decic(odefun,t0,y0,fixed_y0,yp0,fixed_yp0,options) は、AbsTol および RelTol の値を指定するオプション構造体 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^{\prime }-y_2\\ 0& =y_1+y_2\end{array}$$

これらの方程式は十分わかりやすいため、変数の矛盾のない初期条件を簡単に読み取ることができます。たとえば、$$y_1=1$$ を固定すると、2 番目の方程式から $$y_2=-1$$ となり、最初の方程式から $$y_1^{\prime }=-1/2$$ となります。$$y_1$$、$$y_1^{\prime }$$、および $$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 = 2×1

     1
    -1

yp0 = 2×1

   -0.5000
         0

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

入力引数

odefun — 解を求める関数

関数ハンドル

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

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

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

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

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

$$\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;

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

例: @myFcn

データ型: function_handle

t0 — 初期時刻

スカラー

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

データ型: single | double

y0 — y 成分の初期推定値

ベクトル

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

データ型: single | double

fixed_y0 — 固定する y 成分

1 と 0 からなるベクトル | []

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

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

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

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

yp0 — y' 成分の初期推定値

ベクトル

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

データ型: single | double

fixed_yp0 — 固定する y' 成分

1 と 0 からなるベクトル | []

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

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

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

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

options — オプション構造体

構造体配列

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

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

データ型: struct

出力引数

y0_new — y0 の矛盾のない初期条件

ベクトル

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

yp0_new — yp0 の矛盾のない初期条件

ベクトル

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

resnrm — 残差のノルム

ベクトル

残差のノルム。ベクトルとして返されます。resnrm は odefun(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 を指定し、odeset の InitialSlope オプションの値として yp_new を指定します。

バージョン履歴

R2006a より前に導入

参考

ode15i | odeget | odeset