ode45

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

次の ODE を解きます。

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

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

解をプロットします。

plot(t,y,'-o')

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

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

関数ファイル 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')

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

次の ODE を解きます。

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

この例の最後にあるローカル関数 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

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 を解くのに便利です。ただし、この手法にはいくつかのトレードオフもあります。

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

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

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: [1x1 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 ... ] (1x60 double)
          y: [2x60 double]
      stats: [1x1 struct]
      idata: [1x1 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')