このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

odeFunction

ODE ソルバー用関数ハンドルへのシンボリック式の変換

構文

f = odeFunction(expr,vars)

f = odeFunction(expr,vars,p1,...,pN)

f = odeFunction(___,Name,Value)

説明

f = odeFunction(expr,vars) は、シンボリック代数式系を MATLAB^® 関数ハンドルに変換します。この関数ハンドルは、ode15i を除く数値 MATLAB ODE ソルバーの入力引数として使用できます。引数 vars は系の状態変数を指定します。

f = odeFunction(expr,vars,p1,...,pN) では、系のシンボリック パラメーターを p1,...,pN に指定します。

f = odeFunction(___,Name,Value) は、1 つ以上の Name,Value 引数のペアによって指定された追加オプションを使用します。

例

ODE ソルバー用関数ハンドルの作成と DAE の求解

シンボリックな微分代数方程式の系を MATLAB ODE ソルバーに適した関数ハンドルに変換します。次に ode15s ソルバーを使用して、この系の解を求めます。

次の 2 階微分代数方程式を作成します。

syms y(t);
eqn = diff(y(t),t,2) == (1-y(t)^2)*diff(y(t),t) - y(t);

reduceDifferentialOrder を使用して、この微分方程式を 2 つの 1 階微分方程式の方程式系として書き換えます。ここで、vars は系の状態変数ベクトルです。新規の変数 Dy(t) は、t に対する y(t) の 1 次導関数を表します。

[eqs,vars] = reduceDifferentialOrder(eqn,y(t))

eqs =
 diff(Dyt(t), t) + y(t) + Dyt(t)*(y(t)^2 - 1)
                       Dyt(t) - diff(y(t), t)
 
vars =
   y(t)
 Dyt(t)

y(t) およびその導関数 Dy(t) の初期条件を 2 および 0 にそれぞれ設定します。

initConditions = [2 0];

系の質量行列 M および方程式 F の右辺を含むベクトルを求めます。

[M,F] = massMatrixForm(eqs,vars)

M =
[  0, 1]
[ -1, 0]
 
F =
 - y(t) - Dyt(t)*(y(t)^2 - 1)
                      -Dyt(t)

M および F は次の形で表されます。$$M\left(t,x\left(t\right)\right)\dot{x}\left(t\right)=F\left(t,x\left(t\right)\right).$$。後の計算を単純化するために、系を次の形に書き換えます。$$\dot{x}\left(t\right)=f\left(t,x\left(t\right)\right)$$。

f = M\F

f =
                          Dyt(t)
 - Dyt(t)*y(t)^2 - y(t) + Dyt(t)

odeFunction を使用して、f を MATLAB 関数ハンドルに変換します。結果の関数ハンドルは MATLAB ODE ソルバー ode15s の入力になります。

odefun = odeFunction(f,vars);
ode15s(odefun, [0 10], initConditions)

シンボリック パラメーターが含まれている系の関数ハンドル

状態変数とシンボリック パラメーターが含まれているシンボリック微分方程式の系を MATLAB ODE ソルバーに適した関数ハンドルに変換します。

微分代数方程式系を作成します。ここで、シンボリック関数 x1(t) および x2(t) は系の状態変数を表しています。また、系には定数シンボリック パラメーター a、b、およびパラメーター関数 r(t) が含まれます。これらのパラメーターは状態変数を表していません。方程式と状態変数を 2 つのシンボリック ベクトル (方程式をシンボリック方程式のベクトル、変数をシンボリック関数呼び出しのベクトル) として指定します。

syms x1(t) x2(t) a b r(t)
eqs = [diff(x1(t),t) == a*x1(t) + b*x2(t)^2,...
       x1(t)^2 + x2(t)^2 == r(t)^2];
vars = [x1(t) x2(t)];

質量行列 M と、この系の右辺から成るベクトル F を求めます。M および F は次の形で表されます。$$M\left(t,x\left(t\right)\right)\dot{x}\left(t\right)=F\left(t,x\left(t\right)\right).$$

[M,F] = massMatrixForm(eqs,vars)

M =
[ 1, 0]
[ 0, 0]
 
F =
        b*x2(t)^2 + a*x1(t)
 r(t)^2 - x1(t)^2 - x2(t)^2

odeFunction を使用して、M および F から MATLAB 関数ハンドルを生成します。関数ハンドル F にはシンボリック パラメーターが含まれています。

M = odeFunction(M,vars)
F = odeFunction(F,vars,a,b,r(t))

M = 
  function_handle with value:
    @(t,in2)reshape([1.0,0.0,0.0,0.0],[2,2])

F = 
  function_handle with value:
    @(t,in2,param1,param2,param3)[param1.*in2(1,:)+...
    param2.*in2(2,:).^2;param3.^2-in2(1,:).^2-in2(2,:).^2]

パラメーター値を指定します。

a = -0.6;
b = -0.1;
r = @(t) cos(t)/(1+t^2);

簡約された関数ハンドル F を作成します。

F = @(t,Y) F(t,Y,a,b,r(t));

DAE 系に対し整合性のある初期条件を指定します。

t0 = 0;
y0 = [-r(t0)*sin(0.1); r(t0)*cos(0.1)];
yp0 = [a*y0(1) + b*y0(2)^2; 1.234];

系の質量行列 M および、導関数の初期条件のベクトル yp0 を含むオプション セットを作成します。

opt = odeset('mass',M,'InitialSlope',yp0);

次に、ode15s を使用して方程式系の解を求めます。

ode15s(F, [t0, 1], y0, opt)

関数ハンドルとコメントのファイルへの書き込み

File オプションを使用して、生成した関数ハンドルをファイルに書き込みます。ファイルへの書き込みの際に、odeFunction は、t0、t1 といった名前の中間変数を使用してコードを最適化します。Comments オプションを指定してファイルにコメントを含めます。

微分方程式系の定義質量行列 M とその右辺 F を求めます。

syms x(t) y(t)
eqs = [diff(x(t),t)+2*diff(y(t),t) == 0.1*y(t), ...
       x(t)-y(t) == cos(t)-0.2*t*sin(x(t))];
vars = [x(t) y(t)];
[M,F] = massMatrixForm(eqs,vars);

M と F の MATLAB コードをファイル myfileM と myfileF に書き込みます。odeFunction は既存のファイルを上書きします。ファイル内にコメント Version: 1.1 を含めます。出力ファイルは開くことと編集することができます。

M = odeFunction(M,vars,'File','myfileM','Comments','Version: 1.1');

function expr = myfileM(t,in2)
%MYFILEM
%    EXPR = MYFILEM(T,IN2)

%    This function was generated by the Symbolic Math Toolbox version 7.3.
%    01-Jan-2017 00:00:00

%Version: 1.1
expr = reshape([1.0,0.0,2.0,0.0],[2, 2]);

F = odeFunction(F,vars,'File','myfileF','Comments','Version: 1.1');

function expr = myfileF(t,in2)
%MYFILEF
%    EXPR = MYFILEF(T,IN2)

%    This function was generated by the Symbolic Math Toolbox version 7.3.
%    01-Jan-2017 00:00:00

%Version: 1.1
x = in2(1,:);
y = in2(2,:);
expr = [y.*(1.0./1.0e1);-x+y+cos(t)-t.*sin(x).*(1.0./5.0)];

x(t)、y(t) およびそれらの 1 次導関数に矛盾のない初期値を指定します。

xy0 = [2; 1];    % x(t) and y(t)
xyp0 = [0; 0.05*xy0(2)];    % derivatives of x(t) and y(t)

質量行列 M、初期条件 xyp0、および数値検索の数値許容誤差を服務オプション セットを作成します。

opt = odeset('mass', M, 'RelTol', 10^(-6),...
               'AbsTol', 10^(-6), 'InitialSlope', xyp0);

ode15s を使用して方程式系を解きます。

ode15s(F, [0 7], xy0, opt)

スパース行列

スパース シンボリック行列を MATLAB 関数ハンドルに変換するときには、名前と値のペアの引数 'Sparse',true を使用します。

微分代数方程式系を作成します。ここで、シンボリック関数 x1(t) および x2(t) は系の状態変数を表しています。方程式と状態変数を 2 つのシンボリック ベクトル (方程式をシンボリック方程式のベクトル、変数をシンボリック関数呼び出しのベクトル) として指定します。

syms x1(t) x2(t)

a = -0.6;
b = -0.1;
r = @(t) cos(t)/(1 + t^2);

eqs = [diff(x1(t),t) == a*x1(t) + b*x2(t)^2,...
       x1(t)^2 + x2(t)^2 == r(t)^2];
vars = [x1(t) x2(t)];

質量行列 M と、この系の右辺から成るベクトル F を求めます。M および F は次の形で表されます。$$M\left(t,x\left(t\right)\right)\dot{x}\left(t\right)=F\left(t,x\left(t\right)\right).$$

[M,F] = massMatrixForm(eqs,vars)

M =
[ 1, 0]
[ 0, 0]
 
F =
               - (3*x1(t))/5 - x2(t)^2/10
 cos(t)^2/(t^2 + 1)^2 - x1(t)^2 - x2(t)^2

MATLAB 関数ハンドルを M と F から生成します。質量行列 M の要素のほとんどはゼロのため、M を変換するときには引数 Sparse を使用します。

M = odeFunction(M,vars,'Sparse',true)
F = odeFunction(F,vars)

M = 
  function_handle with value:
    @(t,in2)sparse([1],[1],[1.0],2,2)

F = 
  function_handle with value:
    @(t,in2)[in2(1,:).*(-3.0./5.0)-in2(2,:).^2./1.0e+1;...
    cos(t).^2.*1.0./(t.^2+1.0).^2-in2(1,:).^2-in2(2,:).^2]

DAE 系に対し整合性のある初期条件を指定します。

t0 = 0;
y0 = [-r(t0)*sin(0.1); r(t0)*cos(0.1)];
yp0= [a*y0(1) + b*y0(2)^2; 1.234];

系の質量行列 M および、導関数の初期条件のベクトル yp0 を含むオプション セットを作成します。

opt = odeset('mass',M,'InitialSlope', yp0);

ode15s を使用して方程式系を解きます。

ode15s(F, [t0, 1], y0, opt)

入力引数

expr — 代数式の系

シンボリック式のベクトル

代数式の系。シンボリック式のベクトルとして指定します。

vars — 状態変数

シンボリック関数のベクトル | シンボリック関数呼び出しのベクトル

状態変数。x(t) など、シンボリックな関数または関数呼び出しのベクトルとして指定します。

例: [x(t),y(t)] または [x(t);y(t)]

p1,...,pN — 系のパラメーター

シンボリック変数 | シンボリック関数 | シンボリック関数呼び出し | シンボリック ベクトル | シンボリック行列

系のパラメーター。f(t) などのシンボリック変数、シンボリック関数またはシンボリック関数呼び出しで指定します。また、系のパラメーターをシンボリックな変数、関数または関数呼び出しのベクトルまたは行列として指定することもできます。expr に vars で指定されている変数以外のシンボリック パラメーターが含まれる場合は、それらの追加パラメーターを p1,...,pN として指定しなければなりません。

名前と値の引数

オプションの引数のペアを Name1=Value1,...,NameN=ValueN として指定します。ここで、Name は引数名、Value は対応する値です。名前と値の引数は他の引数の後になければなりませんが、ペアの順序は関係ありません。

R2021a より前では、コンマを使用して名前と値の各ペアを区切り、Name を引用符で囲みます。

例: odeFunction(expr,vars,'File','myfile')

Comments — ファイル ヘッダーに含めるコメント

文字ベクトル | 文字ベクトルの cell 配列 | string ベクトル

ファイル ヘッダーに含めるコメント。文字ベクトル、文字ベクトルの cell 配列、または string ベクトルとして指定します。

File — 生成されたコードを含むファイルへのパス

文字ベクトル

生成されたコードを含むファイルへのパス。文字ベクトルとして指定します。生成されたファイルは double 型の引数を受け入れます。また、Symbolic Math Toolbox™ なしで使用できます。値が空の場合、odeFunction は無名関数を生成します。文字ベクトルが .m で終わっていない場合は、関数は .m を追加します。

既定では、引数 File をもつ odeFunction は、最適化されたコードを含むファイルを生成します。"最適化された" とは、コードを単純化または高速化するために中間変数が自動的に生成されていることを意味します。MATLAB によって生成される中間変数は、小文字の t に自動生成される番号を付けたもの (t32 など) です。コードの最適化を無効にするには、引数 Optimize を使用します。

Optimize — 関数ファイルに書き込まれたコードの最適化を回避するフラグ

true (既定値) | false

関数ファイルに書き込まれたコードの最適化を回避するフラグ。false または true で指定します。

既定では、引数 File をもつ odeFunction は、最適化されたコードを含むファイルを生成します。"最適化された" とは、コードを単純化または高速化するために中間変数が自動的に生成されていることを意味します。MATLAB によって生成される中間変数は、小文字の t に自動生成される番号を付けたもの (t32 など) です。

引数 File を持たない (または空の文字ベクトルによって指定されたファイル パスを持つ) odeFunction は、関数ハンドルを作成します。この場合、コードは最適化されません。Optimize を true に設定して強制的にコードの最適化を行おうとすると、odeFunction はエラーをスローします。

Sparse — スパース行列と密行列の生成を切り替えるフラグ

false (既定値) | true

スパース行列と密行列の生成を切り替えるフラグ。true または false で指定します。'Sparse',true を指定すると、生成された関数では、スパース数値行列でシンボリック行列を表します。ゼロの要素を多く含むシンボリック行列を変換する場合、'Sparse',true を指定します。多くの場合、スパース行列での演算は、密行列での同じ演算よりも効率的です。スパース行列を参照してください。

出力引数

f — 関数ハンドルは、ode15i を除く数値 MATLAB ODE ソルバーの入力引数です。

MATLAB 関数ハンドル

ode15i を除くすべての数値 MATLAB ODE ソルバー入力引数として機能する関数ハンドル。MATLAB 関数ハンドルとして返されます。

odeFunction は、ode45、ode15s、ode23t などの ODE ソルバーに適した関数ハンドルを返します。この関数ハンドルを受け付けない唯一の ODE ソルバーは、完全な陰的微分方程式向けのソルバーである ode15i です。方程式系を ode15i に適した関数ハンドルに変換するには、daeFunction を使用します。

バージョン履歴

R2015a で導入

参考

findDecoupledBlocks | daeFunction | decic | incidenceMatrix | isLowIndexDAE | massMatrixForm | matlabFunction | ode15i | ode15s | ode45 | ode23t | reduceDAEIndex | reduceDAEToODE | reduceDifferentialOrder | reduceRedundancies

トピック

• 質量行列ソルバーを使用した DAE の求解