MATLAB Function Reference | ![]() ![]() |
表示
[T,Y] = solver(odefun,tspan,y0) [T,Y] = solver(odefun,tspan,y0,options) [T,Y] = solver(odefun,tspan,y0,options,p1,p2...) [T,Y,TE,YE,IE] = solver(odefun,tspan,y0,options)
ここで、solver
は、ode45
, ode23
, ode113
, ode15s
, ode23s
, ode23t
, ode23tb
のいずれかです。
引数
|
微分方程式の右辺を計算する関数。すべてのソルバは、![]() ![]() ode23s のみは、定質量行列をもつ方程式を解くことができます。ode15s と ode23t は、正則でない質量行列をもつ方程式、すなわち、微分代数方程式を解くことができます。 |
tspan |
積分区間[t0 tfinal] を指定するベクトル。特定の時間(単純増加または単純減少)で解を得るためには、tspan = [t0,t1,...,tfinal] を使います。
|
y0 |
初期条件ベクトル |
|
関数odeset を使って、作成されるオプションの積分引数。odeset を参照してください。 |
p1,p2... |
odefun に渡されるオプションのパラメータ |
詳細
[T,Y] =
は、solver
(odefun,tspan,y0)
tspan = [t0 tf]
のとき、初期状態y0
で、t0
からtfinal
まで微分方程式システム を積分します。スカラ
t
、列ベクトル y
に対して、関数 f = odefun(t,y)
は、 に対応する列ベクトル
f
を出力します。解の配列 Y
の各行は、列ベクトル T
に出力される時間に対応します。時間 t0
, t1,...,tf
(単調増加、または、単調減少)での解を求めるためには、tspan = [t0,t1,...,tf]
を使ってください。
[T,Y] =
は、solver
(odefun,tspan,y0,options)
options
は、関数odeset
で作成された引数です(詳細は、 odeset
を参照してください)。一般的に使用されるプロパティは、スカラの相対許容誤差RelTol
(デフォルトでは、1e-3
)と、絶対許容誤差AbsTol
(デフォルトでは、1e-6
)のベクトルを含みます。
[T,Y] =
は、M-ファイルsolver
(odefun,tspan,y0,options,p1,p2...)
F
がコールされたときに、付加パラメータp1,p2,...
を関数 odefun
に渡して、上記のように解きます。オプションを指定しない場合は、options = []
を使ってください。
[T,Y,TE,YE,IE] =
は、イベント関数をコールして、関数 solver
(odefun,tspan,y0,options)
が、ゼロになる位置を見つけながら、上記のように解きます。各イベント関数に対して、ユーザは、積分が、ゼロで終了しているか、そして、ゼロクロスが生じる方向を決めるかを指定することができます。これは、
Events
プロパティを @EVENTS
に設定し、関数 [value
,isterminal
,direction
] = EVENTS
(t
,y
) を作成することで、処理されます。i
番目のイベント関数に対して、
value(i)
は、関数値です。isterminal(i)
= 1
に、その他の場合 0 になります。direction(i) = 0
、イベント関数が増加する部分でのみゼロになる場合 +1
、減少する部分でのみゼロになる場合 -1
になります。TE
, YE
, IE
の中の対応する要素は、各々、イベントが発生した時刻、イベント時刻での解、ゼロになるイベント関数のインデックス i
を出力します。
OutputFcn
プロパティの値として出力関数を設定する場合、ソルバは、各時間ステップの後で、計算された解と共にコールします。4つの出力関数 odeplot
, odephas2
, odephas3
, odeprint
が用意されています。ソルバを出力引数を設定しないでコールすると、デフォルトの odeplot
がコールされて、解が計算されながら、プロットされます。odephas2
と odephas3
は、2次元、3次元の位相平面プロットをそれぞれ作成します。odeprint
は、スクリーン上に解の成分を表示します。デフォルトでは、解のすべての成分は、出力関数に渡され、OutputSel
プロパティの値として、インデックスのベクトルを設定することにより、指定した成分のみを渡すこともできます。たとえば、出力引数を設定しないで、ソルバをコールし、OutputSel
の値を [1,3]
に設定することにより、ソルバは、解を計算しながら、解の成分1と3をプロットします。
ステッフなソルバ ode15s
, ode23s
, ode23t
, ode23tb
に対して、ヤコビアン行列 は、信頼性と効率に敏感に影響します。
FJAC(T,Y)
が、ヤコビアン の場合。
odeset
を使って、Jacobian
に @FJAC
を設定します。または、ヤコビアンが定数の場合は、行列 を設定します。
Jacobian
プロパティが設定されていない(デフォルト)場合、 は、有限差分で近似されます。ODE 関数において、
odefun
(T
,[Y1,Y2 ...
]) が、[odefun
(T
,Y1
),odefun
(T
,Y2
) ...
] を出力するようにコード化されている場合、Vectorized
プロパティを 'on
' に設定します。 がスパース行列の場合、
JPattern
プロパティを のスパース性パターン、すなわち、
の
i
番目の要素が、 の
j
番目の要素に依存する場合、S(i,j)
= 1 を、その他の場合、S(i,j)
= 0 を示すパターンを設定します。
ODE suite のソルバは、時間と状態に依存する質量行列 をもつ問題
を解くことができます(ソルバ
ode23s
は、定質量行列をもつ方程式のみを解くことができます)。問題が質量行列をもつ場合、質量行列の値を出力する関数 を作成し、
odeset
を使って、Mass
プロパティを @MASS
に設定します。質量行列が定数の場合、行列が、Mass
プロパティの値として使われます。状態依存の質量行列をもつ問題は、より難しいものになります。
MASS
が、一入力引数 t
と共にコールされる場合、MStateDependence
プロパティを 'none
' に設定します。MStateDependence
を 'weak
' (デフォルト) に設定し、他の場合は、'strong
' に設定してください。どちらの場合でも、関数
MASS
は、2つの引数 (t
,y
) と共にコールされます。多くの微分方程式が存在する場合、スパース性を利用することは重要です。
JPattern
プロパティを使って、Jacobian
プロパティを使って、スパース k
に対して、i,k
) 成分が、j
に依存する場合、S(i,j) = 1
で、その他の場合、S(i,j) = 0
であるスパース行列 S
を MvPattern
に設定します。 質量行列 が正則でない場合、
は、微分代数方程式になります。
が、無矛盾のときにのみ、
であるベクトルが存在すれば、DAEs は、解をもちます。ソルバ
ode15s
と ode23t
は、y0
が、無矛盾である状態に十分に近い場合、インデックス1の DAEs を解くことができます。質量行列が存在する場合、関数 odeset
を使って、MassSingular
プロパティを 'yes'
, 'no'
, 'maybe'
のいずれかに設定できます。デフォルトは、'maybe'
で、問題が、DAE であるか否かのテストを行います。InitialSlope
プロパティの初期値として、yp0
を与えることもできます。デフォルトは、ゼロベクトルです。問題が DAE で、y0
と yp0
が、無矛盾の場合、ソルバは、それらを推定値として取り扱い、推定値に近い無矛盾値を計算しようとし、問題を解き続けます。DAEs を解く場合、 が、対角行列である問題に、定式化できる利点があります。
ODE ソルバの中で使われるアルゴリズムは、精度の次数 [6] とシステムのタイプ(ステッフ、または、ノンステッフ)に従って、解法を設計します。詳細は、アルゴリズム を参照してください。
オプション
種々のソルバは、オプションリストの中で、種々のパラメータを受け入れます。詳細については、odeset
と MATLAB ドキュメントの数学的な解析の微分方程式の ODE ソルバ性能の改良 を参照してください。
例題
例題 1. この例題は、ノンスティッフシステムの例で、外力のない剛体の挙動についての方程式システムです。
このシステムをシミュレーションするために、つぎの方程式を含むファンクションM-ファイル 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);
plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')
例題 2.振動に関するvan der Pol 式によるスティッフなシステムの例です。このリミットサイクルには、解の成分がゆっくりと変化し、問題がスティッフな領域とスティッフでない非常に激しい変化のある領域が交互にあります。
このシステムをシミュレーションするために、つぎの方程式を含むファンクションM-ファイル 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-3、1e-6
)を使い、時間 0 での初期条件ベクトル [2 0] を用いて、時間区間[0 3000] で解きます。
[T,Y] = ode15s(@vdp1000,[0 3000],[2 0]);
出力行列 Y
の最初の列とT
をプロットして、解を示します。
plot(T,Y(:,1),'-o')
:
アルゴリズム
ode45
は、Dormand-Prince の Runge-Kutta (4,5)の明示的な式に基づいています。これは、y(t
n
)
の計算における1ステップソルバであり、1つ前の時間での解 y(t
n-1
)
のみを必要とします。一般に、ode45
はほとんどの問題に対して、最初の試みとして行うのに最も適している関数です[3]。
ode23
は、Bogacki と Shampineによる明示的な Runge-Kutta (2,3) を実現したものです。粗い許容値やわずかなスティッフ性のあるときには、ode45
よりも効率的です。ode45
と同様に、ode23
は1ステップソルバです[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式に基づいています。これは、1ステップソルバなので、粗い許容値では、ode15s
よりも効率的かもしれません。これは、ode15s
では、効率的でないスティッフな問題を解くことができます[9]。
ode23t
は、"フリー"な内挿子を使って、台形則を実行します。問題が、わずかなステッフ性をもち、数値的な減衰が存在しない場合、このソルバを使ってください。 ode23t
は、DAEs を解くことができます[10]。
ode23tb
は、TR-BDF2の改良版で、最初のステージで、陰的なRunge-Kutta式を使った台形則で、2番目のステップで、2次の後退微分方程式を使います。過程で、同じ繰り返し行列が、2つのステージで計算に使われます。ode23s
のように、このソルバは、粗い許容誤差においては、ode15s
より、効率的かも知れません[8], [1] 。
参考
@
(function_handle
), odeset
, odeget
参考文献
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.
Bogacki, P. and L. F. Shampine, "A 3(2) pair of Runge-Kutta formulas," Appl. Math. Letters, Vol. 2, 1989, pp 1-9.
Dormand, J. R. and P. J. Prince, "A family of embedded Runge-Kutta formulae," J. Comp. Appl. Math., Vol. 6, 1980, pp 19-26.
Forsythe, G. , M. Malcolm, and C. Moler, Computer Methods for Mathematical Computations, Prentice-Hall, New Jersey, 1977.
Kahaner, D. , C. Moler, and S. Nash, Numerical Methods and Software, Prentice-Hall, New Jersey, 1989.
Shampine, L. F. , Numerical Solution of Ordinary Differential Equations, Chapman & Hall, New York, 1994.
Shampine, L. F. and M. K. Gordon, Computer Solution of Ordinary Differential Equations: the Initial Value Problem, W. H. Freeman, San Francisco, 1975.
Shampine, L. F. and M. E. Hosea, "Analysis and Implementation of TR-BDF2," Applied Numerical Mathematics 20, 1996.
Shampine, L. F. and M. W. Reichelt, "The MATLAB ODE Suite," SIAM Journal on Scientific Computing, Vol. 18, 1997, pp 1-22.
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.
![]() | nzmax | odefile | ![]() |