ドキュメンテーション センター

  • 評価版
  • 製品アップデート

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

ode45

ノンスティッフ微分方程式の求解: 中次数法

構文

[T,Y] = solver(odefun,tspan,y0)
[T,Y] = solver(odefun,tspan,y0,options)
[T,Y,TE,YE,IE] = solver(odefun,tspan,y0,options)
sol = solver(odefun,[t0 tf],y0...)

このページでは、ソルバー関数の概要を説明します。ode23ode45ode113ode15sode23sode23t および ode23tb。プレースホルダー solver を任意の関数名に置き換えることで、以下のいずれのソルバーでも呼び出すことができます。

引数

次の表でソルバーの入力引数を説明します。

odefun

微分方程式の右辺を計算する関数ハンドル。すべてのソルバーは、y′ = f(t,y) の形式の方程式または質量行列 M(t,y)y′ = f(t,y) を含む問題を解くことができます。関数 ode23s ソルバーは、定質量行列をもつ方程式のみを解くことができます。関数 ode15s と関数 ode23t は、正則でない質量行列をもつ方程式、つまり、微分代数方程式 (DAE) を解くことができます。

tspan

積分区間 [t0,tf] を指定するベクトル。ソルバーは、tspan(1) で初期条件を設定し、tspan(1) から tspan(end) まで積分します。特定の時間 (単調増加または単調減少) での解を得るには、tspan = [t0,t1,...,tf] を使用します。

2 要素 [t0 tf] をもつ tspan ベクトルを指定すると、ソルバーは、積分ステップごとに計算された解を返します。3 要素以上の tspan ベクトルを指定すると、ソルバーは与えられた時点で計算された解を返します。時点は、昇順または降順のいずれかでなければなりません。

 

3 要素以上の tspan を指定すると、ソルバーが区間 tspan(1) から tspan(end) に進むときに使用する内部タイム ステップは影響を受けません。ODE 関連のすべてのソルバーは、基本の定式を連続的に拡張する手法で出力値を得ます。しかし、ソルバーは必ずしも tspan で指定された時点を正確にたどりません。指定された時点で出力される解は、内部時点で計算された解と同じ次数の精度になります。

3 要素以上の tspan を指定した場合、計算の効率性に与える影響は小さいですが、大規模なシステムでは、メモリ管理に影響を与えます。

y0

初期条件のベクトル。

options

既定の積分プロパティを変更するオプション パラメーターの構造体。これは、4 番目の入力引数です。

[t,y] = solver(odefun,tspan,y0,options)

関数 odeset を使用して、オプションを作成することができます。詳細は、odeset を参照してください。

次の表に、ソルバーの出力引数を示します。

T

時点の列ベクトル。

Y

解の配列。Y の各行は、T の対応する行に返される時点での解です。

TE

イベントが発生した時点。

YE

イベント時点での解。

IE

ゼロになるイベント関数のインデックス i

sol

解を計算する構造体。

説明

tspan = [t0 tf] における [T,Y] = solver(odefun,tspan,y0) は、微分方程式 y′ = f(t,y) により初期条件 y0 で時間 t0 から tf の範囲を積分します。最初の引数 odefun は関数ハンドルです。関数 f = odefun(t,y) (t はスカラー、y は列ベクトル) は、f(t,y) に対応する列ベクトル f を返さなければなりません。解の配列 Y の各行は、列ベクトル T に返される時間に対応します。特定の時間 t0t1,...,tf (単調増加または単調減少) での解を得るには、tspan = [t0,t1,...,tf] を使用します。

「関数のパラメーター化」では、必要に応じて関数 fun にパラメーターを追加する方法について説明しています。

[T,Y] = solver(odefun,tspan,y0,options) は、options に指定された「プロパティ値」で既定の積分パラメーターを置き換えて、上記のように解きます。options は、関数 odeset で作成された引数です。一般的に使用されるプロパティは、スカラーの相対許容誤差 RelTol (既定の設定は 1e-3) と、絶対許容誤差 AbsTol (既定の設定は、すべての成分で 1e-6) のベクトルを含みます。解のある要素が非負でなければならない場合、関数 odeset を使用して、これらの要素のインデックスに、NonNegative プロパティを設定します。詳細は、odeset を参照してください。

[T,Y,TE,YE,IE] = solver(odefun,tspan,y0,options) は、呼び出されたイベント関数である (t,y) の関数がゼロになる位置を見つけながら、上記のように解きます。各イベント関数に、積分がゼロで終了しているかと、ゼロクロッシングが生じる方向を決めるかを指定することができます。これは、'Events' プロパティを関数、つまり events または @events に設定し、関数 [value,isterminal,direction] = events(t,y) を作成することで処理されます。eventsi 番目のイベント関数の各出力引数は以下のようになります。

  • value(i) は、関数値です。

  • 積分が、このイベント関数の零点で終了する場合は isterminal(i) = 1 に、その他の場合は 0 になります。

  • すべての零点が計算可能な場合 (既定の設定)、direction(i) = 0、イベント関数が増加する部分でのみゼロになる場合は +1、減少する部分でのみゼロになる場合は -1 になります。

TEYEIE の各要素はそれぞれ、イベントが発生した時点、イベント時点での解、ゼロになるイベント関数のインデックス i を返します。

sol = solver(odefun,[t0 tf],y0...) は、区間 [t0,tf] 上の任意の点で解を計算する関数 deval で使用できる構造体を返します。関数 odefun は、関数ハンドルとして渡す必要があります。構造体 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 オプションに指定された関数により返されるベクトルのインデックス。値は、どのイベントをソルバーが検出したかを示します。

OutputFcn プロパティの値として出力関数を設定する場合、ソルバーは各タイム ステップの後で、計算された解と共に呼び出します。4 つの出力関数 odeplotodephas2odephas3odeprint が用意されています。ソルバーを出力引数を設定しないで呼び出すと、既定の関数 odeplot が呼び出され、解が計算されながらプロットされます。関数 odephas2 と関数 odephas3 はそれぞれ 2次元、3次元の位相平面プロットを作成します。関数 odeprint は、画面上に解の成分を表示します。既定の設定では、ODE ソルバーは、すべての解の成分を出力関数に渡します。OutputSel プロパティの値に、インデックスのベクトルを設定することにより、指定した成分のみを渡すこともできます。たとえば、出力引数を設定せずに、ソルバーを呼び出し、OutputSel の値を [1,3] に設定すると、ソルバーは解を計算しながら解の成分 1 と 3 をプロットします。

スティッフなソルバーの関数 ode15sode23sode23tode23tb の場合、ヤコビ行列 ∂f/∂y は信頼性と効率性に敏感に影響します。odeset を使用して Jacobian を、FJAC(T,Y) がヤコビ ∂f/∂y を返す場合は @FJAC に設定し、ヤコビが定数の場合は行列 ∂f/∂y に設定します。Jacobian プロパティが設定されていない (既定の設定) 場合、∂f/∂y は有限差分で近似されます。ODE 関数は、odefun(T,[Y1,Y2 ...]) が [odefun(T,Y1),odefun(T,Y2) ...] を返すようにコード化されている場合、Vectorized プロパティを 'on' に設定します。∂f/∂y がスパース行列の場合、JPattern プロパティを ∂f/∂y のスパース パターン、つまり f(t,y) の i 番目の要素が、y の j 番目の要素に依存する場合は、スパース行列 S に対して S(i,j) = 1 を、その他の場合は 0 を示すパターンを設定します。

ODE スイートのソルバーは、時間と状態に依存する質量行列 M をもつ M(t,y)y′ = f(t,y) 形式の問題を解くことができます。ソルバーの関数 ode23s は、定質量行列をもつ方程式のみを解くことができます。問題が質量行列をもつ場合、質量行列の値を返す関数 M = MASS(t,y) を作成し、関数 odeset を使用して、Mass プロパティを @MASS に設定します。質量行列が定数の場合、行列が Mass プロパティの値として使用されます。状態依存の質量行列をもつ問題は、より難しいものになります。

  • 質量行列が状態変数 y に依存せず、関数 MASS が 1 入力引数 t と共に呼び出される場合、MStateDependence プロパティを 'none' に設定します。

  • 質量行列の y への依存性が弱い場合は MStateDependence を 'weak' (既定の設定) に設定し、他の場合は 'strong' に設定してください。どちらの場合でも、関数 MASS は、2 つの引数 (t,y) と共に呼び出されます。

微分方程式が多数存在する場合、以下のようなスパース性の利用が重要になります。

  • スパース M(t,y) を返します。

  • JPattern プロパティを使用して、∂f/∂y のスパース パターンを与えるか、Jacobian プロパティを使用して、スパース ∂f/∂y を与えます。

  • 非常に強い状態依存の M(t,y) の場合、任意の k に対して、M(t,y) の (i,k) 成分が y の成分 j に依存する場合は MvPattern を、S(i,j) = 1 をもつスパース行列 S に設定し、その他の場合は 0 に設定します。

質量行列 M が正則でない場合、M(t,y)y′ = f(t,y) は微分代数方程式になります。DAE が解をもつのは、y0 が定数のときのみです。すなわち、M(t0,y0)yp0 = f(t0,y0) であるようなベクトル yp0 がある場合です。ソルバーの関数 ode15sode23t は、y0 が矛盾がない状態に十分に近い場合、インデックス 1 の DAE を解くことができます。質量行列が存在する場合、関数 odeset を使用して、MassSingular プロパティを 'yes''no''maybe' のいずれかに設定できます。既定値は、'maybe' で、問題が、DAE であるかどうかのテストを行います。InitialSlope プロパティの値として、yp0 を与えることもできます。既定の設定はゼロ ベクトルです。問題が DAE で、y0yp0 が矛盾する場合、ソルバーはそれらを推定値として取り扱い、推定値に近い矛盾のない値を計算しようとし、問題を解き続けます。DAE を解く場合、M が対角行列 (半明示的 DAE) をもつように問題を定式化できる利点があります。

ソルバー

問題のタイプ

精度の次数

使用時

ode45

ノンスティッフ

多くの場合使用。ユーザーが最初に使用するソルバー。

ode23

ノンスティッフ

粗い誤差許容量を使用する場合、または中程度にスティッフな問題を解く場合。

ode113

ノンスティッフ

低から高

厳密な誤差許容量を使用する場合、または計算量の多い ODE ファイルを解く場合。

ode15s

スティッフ

低から中

関数 ode45 がスティッフな問題のために遅い場合。

ode23s

スティッフ

粗い誤差許容量を使用して、スティッフな問題を解く場合や質量行列が定数の場合。

ode23t

少しスティッフ

問題が、中程度にスティッフで、数値的な減衰なしの解を必要とする場合。

ode23tb

スティッフ

粗い誤差許容量を使用して、スティッフ システムを解く場合。

ODE ソルバーで使用されるアルゴリズムは、精度の次数 [6] とシステムのタイプ (スティッフ、またはノンスティッフ) に従って解法を設計します。詳細は、「アルゴリズム」を参照してください。

オプション

さまざまなソルバーが、オプション引数にさまざまなパラメーターを受け入れます。詳細は、『MATLAB® 数学』ドキュメンテーションの「odeset」および「積分器のオプション」を参照してください。

パラメーター

ode45

ode23

ode113

ode15s

ode23s

ode23t

ode23tb

RelTol, AbsTol, NormControl

OutputFcn, OutputSel, Refine, Stats

NonNegative

√ *

√ *

√ *

Events

MaxStep, InitialStep

Jacobian, JPattern, Vectorized

Mass
MStateDependence
MvPattern
MassSingular






















InitialSlope

MaxOrder, BDF

    メモ:    問題に質量行列がない場合のみ、関数 ode15sode23tode23tbNonNegative パラメーターを使用できます。

例 1

ノンスティッフ システムの例は、外力のない剛体の動作に関する方程式です。

このシステムをシミュレーションするために、上記の方程式を含む関数 rigid を作成します。

function dy = rigid(t,y)
dy = zeros(3,1);    % a column vector
dy(1) = y(2) * y(3);
dy(2) = -y(1) * y(3);
dy(3) = -0.51 * y(1) * y(2);

この例では、コマンド odeset で許容誤差を変更し、時間 0 での初期条件ベクトル [0 1 1] を使って、[0 12] の区間で解きます。

options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);
[T,Y] = ode45(@rigid,[0 12],[0 1 1],options);

出力配列 Y の各列と T をプロットして、解を表示します。

plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')

例 2

スティッフ システムの例は、緩和振動に関する van der Pol 式によるものです。このリミット サイクルには、解の成分がゆっくりと変化し、問題がスティッフな領域と、スティッフでない非常に激しい変化のある領域が交互にあります。

このシステムをシミュレーションするために、上記の方程式を含む関数 vdp1000 を作成します。

function dy = vdp1000(t,y)
dy = zeros(2,1);    % a column vector
dy(1) = y(2);
dy(2) = 1000*(1 - y(1)^2)*y(2) - y(1);

この問題では、既定の相対誤差と絶対誤差 (1e-31e-6) を使用し、時間 0 での初期条件ベクトル [2 0] を使って、時間区間 [0 3000] で解きます。

[T,Y] = ode15s(@vdp1000,[0 3000],[2 0]);

出力行列 Y の最初の列と T をプロットして、解を表示します。

plot(T,Y(:,1),'-o')

例 3

この例では、時間に依存する項をもつ常微分方程式を解きます。

時間依存パラメーター (2 つのベクトルから与えられるデータ点のみによって定義) をもつ、以下の ODE を考えます。

y'(t) + f(t)y(t) = g(t)

初期条件は、y(0) = 0 です。ここで、関数 f(t) は、n1 列のベクトル tff によって定義され、関数 g(t) は、m1 列のベクトル tgg によって定義されます。

まず、時間依存パラメーター f(t)g(t) を以下のように定義します。

ft = linspace(0,5,25); % Generate t for f
f = ft.^2 - ft - 3; % Generate f(t)
gt = linspace(1,6,25); % Generate t for g
g = 3*sin(gt-0.25); % Generate g(t)

上記で指定されたデータセットを内挿する関数を記述し、指定した時間で、時間に依存する項の値を取得します。

function dydt = myode(t,y,ft,f,gt,g)
f = interp1(ft,f,t); % Interpolate the data set (ft,f) at time t
g = interp1(gt,g,t); % Interpolate the data set (gt,g) at time t
dydt = -f.*y + g; % Evalute ODE at time t

MATLAB 関数 ode45 内で、最初の入力引数に時間を指定して、導関数 myode.m を呼び出します。

Tspan = [1 5]; % Solve from t=1 to t=5
IC = 1; % y(t=0) = 1
[T Y] = ode45(@(t,y) myode(t,y,ft,f,gt,g),Tspan,IC); % Solve ODE

y(t) を時間の関数としてプロットします。

plot(T, Y);
title('Plot of y as a function of time');
xlabel('Time'); ylabel('Y(t)');

詳細

すべて展開する

アルゴリズム

関数 ode45 は、Dormand-Princean の陽的 Runge-Kutta (4,5) 公式に基づきます。これは、単段階ソルバーです。つまり、y(tn) の計算では、直前の y(tn-1) における解のみ必要とします。一般に、関数 ode45 は、ほとんどの問題で最初の試みとして利用すべき最良の関数です。 [3]

関数 ode23 もまた、Bogacki と Shampine の陽的 Runge-Kutta (2,3) に基づいています。粗い許容範囲で計算するときやわずかにスティッフ性がある場合、関数 ode45 より効率的です。関数 ode45 と同様に関数 ode23 は、単段階ソルバーです。 [2]

関数 ode113 は、可変次数の Adams-Bashforth-Moulton PECE ソルバーです。厳しい許容範囲が設定された問題や ODE ファイル関数の計算に特に時間のかかるような問題では、関数 ode45 より効率的な場合があります。関数 ode113 は、多段ソルバーです。つまり、現在の点の解を求めるために、それ以前の数点を使います。 [7]

上記のアルゴリズムは、ノンスティッフなシステムを解くためのものです。これらが過度に遅い場合は、以下のスティッフなソルバーのいずれかを使用してください。

関数 ode15s は、数値微分式 (numerical differentiation formulas: NDF) に基づく可変次数ソルバーです。オプションとして、NDF よりも一般的に効率的でない各種の後退差分方程式 (backward differentiation formulas: BDF、Gear 法とも呼ばれる) も使います。関数 ode113 と同じように、関数 ode15s も、多段ソルバーです。問題がスティッフである可能性がある場合や、微分代数問題を解く場合や、関数 ode45 では失敗したり、非常に非効率であった場合は、関数 ode15s を実行してみてください。 [9], [10]

関数 ode23s は、2 次の Rosenbrock の公式を改良したものをベースにしています。関数 ode23s は、単段階ソルバーであるため、粗い許容範囲の場合には関数 ode15s よりも効率的な場合があります。関数 ode15s では効果的でないある種のスティッフな問題を解くことができます。 [9]

関数 ode23t は、"フリー" 内挿を利用して台形法を実行します。問題が適度にスティッフで数値的減衰のない解が必要な場合に、このソルバーを利用してください。関数 ode23t は、DAE を解くことができます。 [10]

関数 ode23tb は、TR-BDF2 の改良版で、第 1 ステージで陰的 Runge-Kutta 法を使った台形法ステップを使用し、第 2 ステージで 2 次の後退差分方程式を使用します。この過程で、同じ反復行列が 2 つの段階での実行に利用されます。関数 ode23s のように、このソルバーは粗い許容誤差において関数 ode15s より効率的な可能性があります。 [8], [1]

参考文献

[1] Bank, R. E., W. C. Coughran, Jr., W. Fichtner, E. Grosse, D. Rose, and R. Smith, “Transient Simulation of Silicon Devices and Circuits,” IEEE Trans. CAD, 4 (1985), pp. 436–451.

[2] Bogacki, P. and L. F. Shampine, “A 3(2) pair of Runge-Kutta formulas,” Appl. Math. Letters, Vol. 2, 1989, pp. 321–325.

[3] Dormand, J. R. and P. J. Prince, “A family of embedded Runge-Kutta formulae,” J. Comp. Appl. Math., Vol. 6, 1980, pp. 19–26.

[4] Forsythe, G. , M. Malcolm, and C. Moler, Computer Methods for Mathematical Computations, Prentice-Hall, New Jersey, 1977.

[5] Kahaner, D. , C. Moler, and S. Nash, Numerical Methods and Software, Prentice-Hall, New Jersey, 1989.

[6] Shampine, L. F. , Numerical Solution of Ordinary Differential Equations, Chapman & Hall, New York, 1994.

[7] Shampine, L. F. and M. K. Gordon, Computer Solution of Ordinary Differential Equations: the Initial Value Problem, W. H. Freeman, San Francisco, 1975.

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

[9] Shampine, L. F. and M. W. Reichelt, “The MATLAB ODE Suite,” SIAM Journal on Scientific Computing, Vol. 18, 1997, pp. 1–22.

[10] 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.

参考

| | | |

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