Main Content

decic

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

説明

[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 も使用します。オプション構造体を作成するには odeset を使用します。

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

すべて折りたたむ

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

0=2y1-y20=y1+y2

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

decic を使用し、値 y1=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

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

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

ty2(y)3-y3(y)2+t(t2+1)y-t2y=0.

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

方程式のコード化

ode15i に適した形式で方程式のコードを作成するには、ty、および y への入力を含み、方程式の残差値を返す関数を記述する必要があります。関数 @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 ソルバーは "矛盾のない" 初期条件を必要とします。つまり、ソルバーに提供される初期条件は以下を満たさなければなりません。

f(t0,y,y)=0.

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

ここでは、初期値 y(t0)=32 を固定し、y(t0)=0 の初期推定から計算を開始して、decic で導関数 y(t0) の矛盾のない初期条件を計算します。

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 を時間区間 [1 10] について解きます。

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

結果のプロット

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

y(t)=t2+12.

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

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

Figure contains an axes object. The axes object contains 2 objects of type line. One or more of the lines displays its values using only markers These objects represent ode15i, exact.

入力引数

すべて折りたたむ

解を求める関数。積分する関数を定義する関数ハンドルとして指定します。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 より前に導入

参考

| |