ドキュメンテーション

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

ode23t

中程度にスティッフな常微分方程式 (ODE) と微分代数方程式 (DAE) の求解 — 台形則

構文

  • [t,y] = ode23t(odefun,tspan,y0)
  • [t,y] = ode23t(odefun,tspan,y0,options)
  • [t,y,te,ye,ie] = ode23t(odefun,tspan,y0,options)
  • sol = ode23t(___)

説明

[t,y] = ode23t(odefun,tspan,y0) は、tspan = [t0 tf] のときに、初期条件 y0 を使用して、微分方程式系 y'=f(t,y)t0 から tf まで積分します。解の配列 y の各行は、列ベクトル t に返される値に対応します。

すべての MATLAB® ODE ソルバーは、y'=f(t,y) の形式の方程式系、あるいは質量行列 M(t,y)y'=f(t,y) を含む問題を解くことができます。ode23s ソルバーは、質量行列が定数である場合にのみ、これを含む問題を解くことができます。ode15s および ode23t は、特異質量行列をもつ方程式、つまり微分代数方程式 (DAE) を解くことができます。odesetMass オプションを使用して質量行列を指定します。

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

[t,y,te,ye,ie] = ode23t(odefun,tspan,y0,options) はさらに、(t,y) の関数 (イベント関数) がゼロになる点を求めます。出力の te はイベント時点、ye はイベント時点における解、ie はトリガーされたイベントのインデックスです。

各関数に対して、ゼロで積分を終了するかどうかと、ゼロクロッシングの方向を考慮するかどうかを指定します。これを行うには、myEventFcn@myEventFcn などの関数に 'Events' プロパティを設定し、対応する関数 [value,isterminal,direction] = myEventFcn(t,y) を作成します。詳細については、「ODE のイベント検出」を参照してください。

sol = ode23t(___) は、区間 [t0 tf] の任意の点で解を計算する関数 deval で使用できる構造体を返します。前述の構文にある任意の入力引数の組み合わせが使用できます。

すべて折りたたむ

単一解要素をもつシンプルな ODE は、ソルバーの呼び出し内に無名関数として指定できます。この無名関数は、2 つの入力 (t,y) を受け入れなければなりません (いずれかの入力が使用されない場合でも)。

次の ODE を解きます。

$$y' = -10t.$$

時間区間 [0,2] および初期条件 y0 = 1 を使用します。

tspan = [0 2];
y0 = 1;
[t,y] = ode23t(@(t,y) -10*t, tspan, y0);

解をプロットします。

plot(t,y,'-o')

スティッフな方程式系の例の 1 つに緩和振動に関するファン デル ポールの方程式があります。このリミット サイクルには、解の成分がゆっくり変化して問題がかなりスティッフな領域と、急激に変化してスティッフでない領域が交互にあります。

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

$$\begin{array}{cl} y_1' &= y_2\\y_2' &= 1000(1-y_1^2)y_2-y_1\end{array}$$

初期条件は、 $y_1(0)=2$ および $y_2(0)=0$ です。MATLAB® に付属の関数 vdp1000 が方程式をエンコードします。

function dydt = vdp1000(t,y)
%VDP1000  Evaluate the van der Pol ODEs for mu = 1000.
%
%   See also ODE15S, ODE23S, ODE23T, ODE23TB.

%   Jacek Kierzenka and Lawrence F. Shampine
%   Copyright 1984-2014 The MathWorks, Inc.

dydt = [y(2); 1000*(1-y(1)^2)*y(2)-y(1)];

既定の相対許容誤差および絶対許容誤差 (それぞれ 1e-3 および 1e-6) によって ode45 を使用し、この系を求解すると非常に時間がかかります。また、解のプロットにも数分かかります。許容誤差内に収めるためのスティッフな領域が原因で、ode45 で積分を完了するには数百万ものタイム ステップが必要です。

次のプロットは ode45 で得られた解ですが、この計算には長い時間がかかります。スティッフな領域を通過するために、膨大な数のタイム ステップが必要であることがわかります。

ode23t ソルバーを使用してスティッフな系を解き、解 y の 1 列目を時点 t に対してプロットします。ode23t ソルバーは ode45 よりもはるかに少ないステップでスティッフな領域を通過します。

[t,y] = ode23t(@vdp1000,[0 3000],[2 0]);
plot(t,y(:,1),'-o')

ode23t は、2 つの入力引数 ty を使用する関数にのみ使用できます。しかし、関数の外部で定義した追加のパラメーターを、関数ハンドルを指定するタイミングで渡すことができます。

次の ODE を解きます。

$$y'' = \frac{A}{B} t y.$$

この方程式を 1 次系として書き直すと次のようになります。

$$\begin{array}{cl} y'_1 &= y_2\\ y'_2 &= \frac{A}{B} t y_1.
\end{array}$$

odefcn.m はこの方程式系を 4 つの入力引数 (tyAB) を受け入れる関数として表します。

function dydt = odefcn(t,y,A,B)
dydt = zeros(2,1);
dydt(1) = y(2);
dydt(2) = (A/B)*t.*y(1);

ode23t を使用して ODE を解きます。事前定義された AB の値を odefcn に渡す関数ハンドルを指定します。

A = 1;
B = 2;
tspan = [0 5];
y0 = [0 0.01];
[t,y] = ode23t(@(t,y) odefcn(t,y,A,B), tspan, y0);

結果をプロットします。

plot(t,y(:,1),'-o',t,y(:,2),'-.')

ode15s ソルバーは、ほとんどのスティッフな問題に適する最初の選択肢です。しかし、特定のタイプの問題に対しては他のスティッフ ソルバーの方が効率的な場合もあります。この例では、4 つすべてのスティッフ ODE ソルバーを使用してスティッフなテスト方程式を解きます。

次のテスト方程式を考えます。

$$y' = -\lambda y.$$

この方程式は $\lambda$ の絶対値が大きくなるにつれてスティッフ性が増大します。 $\lambda = 1 \times 10^9$、初期条件 $y(0)=1$ および時間区間 [0 0.5] を使用します。これらの値を使用すると、問題が十分にスティッフであるため、ode45 および ode23 では方程式の積分が困難です。また、odeset を使用して定数ヤコビ行列 $J = \frac{\partial f}{\partial y} = - \lambda$ を渡し、ソルバー統計の表示をオンにします。

lambda = 1e9;
opts = odeset('Stats','on','Jacobian',-lambda);
tspan = [0 0.5];
y0 = 1;
subplot(2,2,1)
disp('ode15s stats:')
tic, ode15s(@(t,y) -lambda*y, tspan, y0, opts), toc
title('ode15s')

subplot(2,2,2)
disp(' ')
disp('ode23s stats:')
tic, ode23s(@(t,y) -lambda*y, tspan, y0, opts), toc
title('ode23s')

subplot(2,2,3)
disp(' ')
disp('ode23t stats:')
tic, ode23t(@(t,y) -lambda*y, tspan, y0, opts), toc
title('ode23t')

subplot(2,2,4)
disp(' ')
disp('ode23tb stats:')
tic, ode23tb(@(t,y) -lambda*y, tspan, y0, opts), toc
title('ode23tb')
ode15s stats:
104 successful steps
1 failed attempts
212 function evaluations
0 partial derivatives
21 LU decompositions
210 solutions of linear systems
Elapsed time is 1.275721 seconds.
 
ode23s stats:
63 successful steps
0 failed attempts
191 function evaluations
0 partial derivatives
63 LU decompositions
189 solutions of linear systems
Elapsed time is 0.421286 seconds.
 
ode23t stats:
95 successful steps
0 failed attempts
125 function evaluations
0 partial derivatives
28 LU decompositions
123 solutions of linear systems
Elapsed time is 0.616641 seconds.
 
ode23tb stats:
71 successful steps
0 failed attempts
167 function evaluations
0 partial derivatives
23 LU decompositions
236 solutions of linear systems
Elapsed time is 0.498976 seconds.

すべてのスティッフ ソルバーは良好に動作しますが、この特定の問題では ode23s が最小ステップ数、最短時間で積分を完了します。定数ヤコビ行列が指定されているため、どのソルバーも解を求めるために偏導関数を計算する必要はありません。通常、ode23s はステップごとにヤコビ行列を評価するため、ヤコビ行列の指定は ode23s に最も利点があります。

一般のスティッフな問題の場合、スティッフ ソルバーのパフォーマンスは問題の形式と指定されたオプションによって異なります。ヤコビ行列やスパース パターンを指定すると、スティッフな問題に対するソルバーの効率が常に向上します。しかし、各スティッフ ソルバーでのヤコビ行列の扱いは異なるため、向上の程度は大きく異なる可能性があります。実際、方程式系が非常に大きい場合、あるいは何回も解く必要がある場合は、実行時間を最短にするためにさまざまなソルバーのパフォーマンスを調べる価値があります。

この例では、ode23t を使用して電気回路に関するスティッフな微分代数方程式 (DAE) 問題を解く方法を説明します。この問題の原典は次のとおりです: E. Hairer and G. Wanner、Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems、2nd ed.、337 ページ (Berlin: Springer, 1996)。

amp1dae.m にコード記述されている 1 トランジスタ アンプ問題は半陽的な形式に書き換えることができますが、この例では元の形式 $M u' = \phi(u)$ で解きます。この問題に含まれる定数の特異質量行列は、対角行列ではありません。このトランジスタ アンプ回路には 6 つの抵抗、3 つのコンデンサ、1 つのトランジスタが含まれます。

初期の電圧信号は $U_e(t)=0.4 \sin\left( 200 \pi t \right)$、他のパラメーターは定数です。目的はノード 5 を通る出力電圧を求めることです。

初期電圧および動作電圧の値を設定します。

Ue = @(t) 0.4*sin(200*pi*t);
Ub = 6;

キルヒホッフの法則を使用して各ノード (1 ~ 5) を通る電流を等しくすると、5 つの方程式からなる系が得られます。この系の質量行列は次の形式になります。

$$ M = \left( \begin{array}{ccccc} -c_1 & c_1 & 0 & 0 & 0\\ c_1 & -c_1 &
0 & 0 & 0\\ 0 & 0 & -c_2 & 0 & 0\\ 0 & 0 & 0 & -c_3 & c_3\\ 0 & 0 & 0 &
c_3 & -c_3 \end{array} \right),$$

ここで、 $k=1,2,3$ について $c_k = k \times 10^{-6}$ です。

関数 odeset を使用して質量行列を渡します。ソルバーによる DAE 問題の自動検出をテストするために、'MassSingular' オプションを既定値 'maybe' のままにします。

c = 1e-6 * (1:3);
M = zeros(5,5);
M(1,1) = -c(1);
M(1,2) =  c(1);
M(2,1) =  c(1);
M(2,2) = -c(1);
M(3,3) = -c(2);
M(4,4) = -c(3);
M(4,5) =  c(3);
M(5,4) =  c(3);
M(5,5) = -c(3);
opts = odeset('Mass',M);

関数 transampdae.m にはこの例の方程式系が含まれています。この関数を現在のフォルダーに保存して、例を実行します。

function dudt = transampdae(t,u)
% Define voltages and constant parameters
Ue = @(t) 0.4*sin(200*pi*t);
Ub = 6;
R0 = 1000;
R15 = 9000;
alpha = 0.99;
beta = 1e-6;
Uf = 0.026;

f23 = beta*(exp((u(2) - u(3))/Uf) - 1);
dudt = [ -(Ue(t) - u(1))/R0
    -(Ub/R15 - u(2)*2/R15 - (1-alpha)*f23)
    -(f23 - u(3)/R15)
    -((Ub - u(4))/R15 - alpha*f23)
    (u(5)/R15) ];
end

初期条件を設定します。この例では Hairer および Wanner が計算した、矛盾のない初期条件を使用します。

u0(1) = 0;
u0(2) = Ub/2;
u0(3) = Ub/2;
u0(4) = Ub;
u0(5) = 0;

ode23t を使用して、時間区間 [0 0.05] での DAE 系の解を求めます。

tspan = [0 0.05];
[t,u] = ode23t(@transampdae,tspan,u0,opts);

初期電圧 $U_e(t)$ と出力電圧 $U_5(t)$ をプロットします。

plot(t,Ue(t),'o',t,u(:,5),'.')
axis([0 0.05 -3 2]);
legend('Input Voltage U_e(t)','Output Voltage U_5(t)', 'Location', 'NorthWest');
title('One transistor amplifier DAE problem solved by ODE23T');
xlabel('t');

入力引数

すべて折りたたむ

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

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

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

function dydt = odefun(t,y)
dydt = 5*y-3;

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

y'1=y1+2y2y'2=3y1+2y2

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

function dydt = odefun(t,y)
dydt = zeros(2,1);
dydt(1) = y(1)+2*y(2);
dydt(2) = 3*y(1)+2*y(2);

関数 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

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

データ型: single | double

オプション構造体。構造体配列として指定します。options 構造体の作成または変更を行うには、関数 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 オプションに指定された関数により返されるベクトルのインデックス。値は、どのイベントをソルバーが検出したかを示します。

詳細

すべて折りたたむ

アルゴリズム

関数 ode23t は、"フリー" 内挿を利用して台形則を実行します。問題が中程度にスティッフであり、数値的減衰のない解が必要な場合には、ode15s よりもこのソルバーが推奨されます。ode23t は微分代数方程式 (DAE) を解くこともできます [1][2]

参照

[1] Shampine, L. F., M. W. Reichelt, and J.A. Kierzenka, “Solving Index-1 DAEs in MATLAB and Simulink”, SIAM Review, Vol. 41, 1999, pp. 538–552.

[2] Shampine, L. F. and M. E. Hosea, “Analysis and Implementation of TR-BDF2,” Applied Numerical Mathematics 20, 1996.

参考

| | |

R2006a より前に導入

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