ode45

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

構文

[t,y] = ode45(odefun,tspan,y0)

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

[t,y,te,ye,ie] = ode45(odefun,tspan,y0,options)

sol = ode45(___)

説明

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

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

ode45 は汎用的な ODE ソルバーであり、ほとんどの問題に対して最初に試すべきソルバーです。ただし、問題がスティッフである場合や高い精度が要求される場合、他の ODE ソルバーの方が適していることもあります。詳細については、ODE ソルバーの選択を参照してください。

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

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

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

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

例

単一解要素をもつ ODE

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

次の ODE を解きます。

$$y^{\prime }=2t.$$

時間区間 [0 5] および初期条件 y0 = 0 を指定します。

tspan = [0 5];
y0 = 0;
[t,y] = ode45(@(t,y) 2*t, tspan, y0);

解をプロットします。

plot(t,y,'-o')

ノンスティッフ方程式の求解

ファン デル ポールの方程式は次のように 2 次 ODE です。

$$y_1^{\prime \prime }-\mu (1-y_1^2)y_1^{\prime }+y_1=0,$$

ここで、$$\mu >0$$ はスカラー パラメーターです。代入 $$y_1^{\prime }=y_2$$ を行って、この方程式を 1 次 ODE 系として書き換えます。その結果得られる 1 次の ODE は、次のようになります。

$$\begin{array}{cl}{y}_1^{\prime }& =y_2\\ y_2^{\prime }& =\mu (1-y_1^2)y_2-y_1.\end{array}$$

関数ファイル vdp1.m は、$$\mu =1$$ を使用するファン デル ポールの方程式を表します。変数 $$y_1$$ と変数 $$y_2$$ は、2 要素ベクトル dydt のエントリ y(1) および y(2) です。

function dydt = vdp1(t,y)
%VDP1  Evaluate the van der Pol ODEs for mu = 1
%
%   See also ODE113, ODE23, ODE45.

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

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

関数 ode45 を使用し、初期値 [2 0] を指定して、時間区間 [0 20] でこの ODE の解を求めます。その結果、出力として時間の列ベクトル t と解の配列 y が得られます。y の各行は、t の対応する行に返される時刻と対応します。y の 1 列目は $$y_1$$ に対応し、2 列目は $$y_2$$ に対応します。

[t,y] = ode45(@vdp1,[0 20],[2; 0]);

$$y_1$$ および $$y_2$$ の解を t に対してプロットします。

plot(t,y(:,1),'-o',t,y(:,2),'-o')
title('Solution of van der Pol Equation (\mu = 1) with ODE45');
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2')

ODE 関数に追加パラメーターを渡す

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

次の ODE を解きます。

$$y^{\prime \prime }=\frac{A}{B}ty.$$

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

$$\begin{array}{cl}{y}_1^{\prime }& =y_2\\ y_2^{\prime }& =\frac{A}{B}t{y}_1.\end{array}$$

この例の最後にあるローカル関数 odefcn は、この連立方程式を 4 つの入力引数 (t、y、A、B) を受け入れる関数として表します。

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

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

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

結果をプロットします。

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

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

複数の初期条件をもつ ODE を解く

1 つの方程式をもつシンプルな ODE 系の場合、複数の初期条件を含むベクトルとして y0 を指定できます。この手法では、スカラー拡張を使用して独立した連立方程式が作成され (それぞれが各初期値に対応)、ode45 は連立方程式の解を求めて各初期値の結果を生成します。

無名関数を作成し、方程式 $\mathit{f}\left(\mathit{t},\mathit{y}\right)=-2\mathit{y}+2\cos \left(\mathit{t}\right)\sin \left(2\mathit{t}\right)$ を表します。関数は、t および y の 2 つの入力を受け入れなければなりません。

yprime = @(t,y) -2*y + 2*cos(t).*sin(2*t);

$\left[-5,5\right]$ の範囲でさまざまな初期条件のベクトルを作成します。

y0 = -5:5; 

ode45 を使用して時間区間 $\left[0,3\right]$ で、各初期条件について方程式の解を求めます。

tspan = [0 3];
[t,y] = ode45(yprime,tspan,y0);

結果をプロットします。

plot(t,y)
grid on
xlabel('t')
ylabel('y')
title('Solutions of y'' = -2y + 2 cos(t) sin(2t), y(0) = -5,-4,...,4,5','interpreter','latex')

この手法は、複数の初期条件をもつシンプルな ODE を解くのに便利です。ただし、この手法にはいくつかのトレードオフもあります。

• 複数の初期条件をもつ連立方程式を解くことはできません。この手法は、複数の初期条件をもつ 1 つの方程式を解く場合にのみ機能します。

• 各ステップでソルバーが選択するタイム ステップは、最小ステップを取るために必要なシステムの方程式に基づいています。つまり、ソルバーは小さなステップを取って 1 つの初期条件について方程式を満たすことはできますが、他の方程式がそれぞれの解を求める場合、別のステップ サイズを使用します。それにもかかわらず、一般的に、複数の初期条件について同時に解を求める方が、for ループを使用して個別に方程式を解くよりも高速です。

この手法の詳細については、複数の初期条件をもつ ODE 系の求解を参照してください。

時間依存の項をもつ ODE

時間依存のパラメーターをもつ次の ODE について考えます。

初期条件は です。関数 f(t) は時点 ft で評価される n 行 1 列のベクトル f で定義されます。関数 g(t) は時点 gt で評価される m 行 1 列のベクトル g で定義されます。

ベクトル f およびベクトル g を作成します。

ft = linspace(0,5,25);
f = ft.^2 - ft - 3;

gt = linspace(1,6,25);
g = 3*sin(gt-0.25);

f と g を内挿し、指定された時点における時間依存の項の値を取得する関数 myode を作成します。この関数を現在のフォルダーに保存して、例の残りの部分を実行します。

関数 myode は各タイム ステップで ODE を評価するための追加の入力引数を受け入れますが、ode45 は最初の 2 つの入力引数 t と y のみを使用します。

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; % Evaluate ODE at time t

ode45 を使用して、時間区間 [1 5] での方程式の解を求めます。ode45 が myode の最初の 2 つの入力引数のみを使用するように、関数ハンドルを使用して関数を指定します。また、odeset を使用して誤差のしきい値を緩和します。

tspan = [1 5];
ic = 1;
opts = odeset('RelTol',1e-2,'AbsTol',1e-4);
[t,y] = ode45(@(t,y) myode(t,y,ft,f,gt,g), tspan, ic, opts);

解 y を時点 t の関数としてプロットします。

plot(t,y)

解の構造体の評価と拡張

ファン デル ポールの方程式は次のように 2 次 ODE です。

$$y_1^{\prime \prime }-\mu (1-y_1^2)y_1^{\prime }+y_1=0.$$

ode45 を使用して、$$\mu =1$$ のファン デル ポールの方程式を解きます。MATLAB® に付属の関数 vdp1.m が方程式をエンコードします。ソルバーや評価点などの解に関する情報をもつ構造体を返すために、1 つの出力を指定します。

tspan = [0 20];
y0 = [2 0];
sol = ode45(@vdp1,tspan,y0)

sol = struct with fields:
     solver: 'ode45'
    extdata: [1×1 struct]
          x: [0 1.0048e-04 6.0285e-04 0.0031 0.0157 0.0785 0.2844 0.5407 0.8788 1.4032 1.8905 2.3778 2.7795 3.1285 3.4093 3.6657 3.9275 4.2944 4.9013 5.3506 5.7998 6.2075 6.5387 6.7519 6.9652 7.2247 7.5719 8.1226 8.6122 9.1017 9.5054 … ] (1×60 double)
          y: [2×60 double]
      stats: [1×1 struct]
      idata: [1×1 struct]

linspace を使用して、区間 [0 20] に 250 点を生成します。deval を使用して、これらの点で解を評価します。

x = linspace(0,20,250);
y = deval(sol,x);

解の最初の要素をプロットします。

plot(x,y(1,:))

odextend を使用して解を $$t_f=35$$ に拡張し、結果を元のプロットに追加します。

sol_new = odextend(sol,@vdp1,35);
x = linspace(20,35,350);
y = deval(sol_new,x);
hold on
plot(x,y(1,:),'r')

入力引数

odefun — 解を求める関数

関数ハンドル

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

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

たとえば、$$y\text{'}=5y-3$$ を解くには、次の関数を使用します。

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

連立方程式の場合、odefun の出力はベクトルです。ベクトルの各要素は 1 つの方程式の右辺の計算値です。たとえば、次の 2 つを含む連立方程式について考えます。

$$\begin{array}{l}y{\text{'}}_1=y_1+2y_2\\ y{\text{'}}_2=3y_1+2y_2\end{array}$$

各タイム ステップで各方程式の右辺の値を計算する関数は以下のとおりです。

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

関数 odefun に追加パラメーターを指定する方法の詳細については、関数のパラメーター化を参照してください。

例: @myFcn

データ型: function_handle

tspan — 積分区間

ベクトル

積分区間。ベクトルとして指定します。少なくとも、tspan は開始時点と終了時点を指定する 2 要素ベクトル [t0 tf] でなければなりません。t0 と tf の間の特定時点における解を取得するには、[t0,t1,t2,...,tf] の形式の長いベクトルを使用します。tspan の要素は、単調増加または単調減少でなければなりません。

ソルバーは開始時点 tspan(1) で y0 によって指定される初期条件を設定し、tspan(1) から tspan(end) まで積分します。

• tspan に 2 つの要素がある場合 ([t0 tf])、ソルバーは区間内の個々の内部積分ステップで計算した解を返します。

• tspan に 2 つを超える要素がある場合 ([t0,t1,t2,...,tf])、ソルバーは指定された各点で計算した解を返します。ただし、ソルバーは tspan で指定された各点に正確にステップするわけではありません。代わりに、ソルバーは独自の内部ステップを使用して解を計算し、tspan 内の要求された点で解を評価します。指定された点で出力された解の精度の次数は、各内部ステップで計算された解と同じです。

  中間点をいくつか指定しても計算効率にほとんど影響しませんが、大規模な方程式の場合はメモリ管理に影響を及ぼす可能性があります。

ソルバーは tspan の値を使用して InitialStep と MaxStep に適した値を計算します。

• tspan に複数の中間点が含まれる場合 ([t0,t1,t2,...,tf])、指定した点は問題のスケールの目安となり、これはソルバーが使用する InitialStep の値に影響することがあります。そのため、ソルバーで得られる解は、tspan を 2 要素ベクトルとして指定するか、中間点を含むベクトルとして指定するかによって異なる場合があります。

• tspan の最初と最後の値は、最大ステップ サイズ MaxStep の計算に使用されます。そのため、tspan の最初または最後の値を変更すると、ソルバーが別のステップ シーケンスを使用する可能性があり、そのため解が変わることがあります。

例: [1 10]

例: [1 3 5 7 9 10]

データ型: single | double

y0 — 初期条件

ベクトル

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

データ型: single | double

options — オプション構造体

構造体配列

オプション構造体。構造体配列として指定します。options 構造体の作成または変更には、関数 odeset を使用します。各 ODE ソルバーと互換性のあるオプションの一覧については、ODE オプションの概要を参照してください。

例: options = odeset('RelTol',1e-5,'Stats','on','OutputFcn',@odeplot) は、相対許容誤差 1e-5 を指定し、ソルバー統計の表示をオンにして、解の計算時に解をプロットする出力関数 @odeplot を指定します。

データ型: struct

出力引数

t — 評価点

列ベクトル

評価点。列ベクトルとして返されます。

• tspan に 2 つの要素がある場合 ([t0 tf])、t には積分の実行に使用される内部評価点が含まれます。

• tspan に 2 つより多い要素が含まれる場合、t は tspan と同じです。

y — 解

配列

解。配列として返されます。y の各行は、t の対応する行に返される値での解です。

te — イベント時点

列ベクトル

イベント時点。列ベクトルとして返されます。te のイベント時点は ye に返された解に対応し、ie は発生したイベントを指定します。

ye — イベント時点での解

配列

イベント時点での解。配列として返されます。te のイベント時点は ye に返された解に対応し、ie は発生したイベントを指定します。

ie — トリガーされるイベント関数のインデックス

列ベクトル

トリガーされるイベント関数のインデックス。列ベクトルとして返されます。te のイベント時点は ye に返された解に対応し、ie は発生したイベントを指定します。

sol — 評価対象の構造体

構造体配列

評価対象の構造体。構造体配列として返されます。この構造体を関数 deval と共に使用して区間 [t0 tf] の任意の点で解を評価します。構造体配列 sol は、常に以下のフィールドを含みます。

構造体フィールド	説明
sol.x	ソルバーによって選択されたステップの行ベクトル。
sol.y	解。各列の sol.y(:,i) には、時点 sol.x(i) での解が含まれます。
sol.solver	ソルバー名。

さらに、odeset の Events オプションを指定してイベントが検出された場合、sol は以下のフィールドも含みます。

構造体フィールド	説明
sol.xe	イベントが発生した点。sol.xe(end) には、終了イベントの厳密な点 (存在する場合) が含まれます。
sol.ye	sol.xe のイベントに対応する解。
sol.ie	Events オプションに指定された関数により返されるベクトルのインデックス。値は、どのイベントをソルバーが検出したかを示します。

アルゴリズム

関数 ode45 は、Dormand-Princean の陽的 Runge-Kutta (4,5) 公式に基づきます。これは、単一ステップ ソルバーです。つまり、y(tn) の計算では、その直前の時点の解 y(tn-1) のみが必要です [1]、[2]。

参照

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

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

拡張機能

C/C++ コード生成

MATLAB® Coder™ を使用して C および C++ コードを生成します。

使用上の注意および制限:

• odeset のオプション引数はすべて定数でなければなりません。

• コード生成では、オプション構造体の中の定数の質量行列はサポートされません。質量行列は関数として指定してください。

• 少なくとも 2 つの出力引数 T と Y を指定しなければなりません。

• 入力タイプは同種でなければなりません (すべて double またはすべて single)。

• 可変サイズのサポートが有効でなければなりません。コード生成では、tspan が 2 つの要素をもつ場合またはイベント関数を使用する場合に、動的メモリ割り当てが必要です。

バージョン履歴

R2006a より前に導入

参考

ode23 | ode78 | ode89 | ode113 | ode15s | odeset | odeget | deval | odextend

トピック

• ODE ソルバーの選択
• ODE オプションの概要
• ノンスティッフ ODE の求解
• 無名関数
• ODE の一般的な問題のトラブルシューティング