定数の遅延をもつ DDE

この例では、dde23 を使用して定数の遅延を伴う DDE (遅延微分方程式) 系の解を求める方法を説明します。

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

$\begin{array}{l} {y_{1}}^{'} (t) = y_{1} (t - 1) \\ {y_{2}}^{'} (t) = y_{1} (t - 1) + y_{2} (t - 0.2) \\ {y_{3}}^{'} (t) = y_{2} (t) . \end{array}$

$t \leq 0$ の履歴関数は定数 $y_{1} (t) = y_{2} (t) = y_{3} (t) = 1$ です。

方程式の時間遅延は $y$ の項にのみ存在し、遅延自体は定数であるため、方程式は "定数遅延" の連立方程式を形成します。

MATLAB® でこの連立方程式を解くには、定数遅延をもつ連立方程式に適した遅延微分方程式ソルバー dde23 を呼び出す前に、方程式、遅延、および履歴をコード化する必要があります。(ここで行ったように) 必要な関数をファイルの最後にローカル関数として含めることも、あるいは個別の名前付きファイルとして MATLAB パスのディレクトリに保存することもできます。

遅延のコード化

最初に、ベクトルを作成して連立方程式の遅延を定義します。この連立方程式には 2 つの異なる遅延があります。つまり、

最初の要素 $y_{1} (t - 1)$ における遅延 1 と、
2 番目の要素 $y_{2} (t - 0.2)$ における遅延 0.2 です。

dde23 は遅延のベクトル引数を受け入れ、その各要素は 1 つの要素の定数遅延に該当します。

lags = [1 0.2];

方程式のコード化

ここで、方程式をコード化する関数を作成します。この関数にはシグネチャ dydt = ddefun(t,y,Z) がなければなりません。ここでは以下のようになります。

t は時間 (独立変数) です。
y は解 (従属変数) です。
Z(:,j) は遅延 $y (t - τ_{j})$ を近似します。ここで、定数遅延 $τ_{j}$ は lags(j) によって与えられます。

これらの入力はソルバーによって自動的に関数に渡されますが、変数名によって方程式のコード化方法が決まります。この場合、

Z(:,1) $\to y_{1} (t - 1)$
Z(:,2) $\to y_{2} (t - 0.2)$

function dydt = ddefun(t,y,Z)
  ylag1 = Z(:,1);
  ylag2 = Z(:,2);

  dydt = [ylag1(1); 
          ylag1(1)+ylag2(2); 
          y(2)];
end

メモ: 関数はすべて例の最後にローカル関数として含まれます。

解の履歴のコード化

次に、解の履歴を定義する関数を作成します。解の履歴は時間 $t \leq t_{0}$ での解です。

function s = history(t)
  s = ones(3,1);
end

方程式の求解

最後に、積分区間 $[t_{0} t_{f}]$ を定義し、dde23 ソルバーを使用して DDE を解きます。

tspan = [0 5];
sol = dde23(@ddefun, lags, @history, tspan);

解のプロット

解の構造体 sol には sol.x フィールドと sol.y フィールドがあり、ソルバーが受け入れる内部タイムステップと、その時間に対応する解が含まれます (特定の点での解が必要な場合、deval を使用して特定の点での解を評価することができます)。

3 つの解要素を時間に対してプロットします。

plot(sol.x,sol.y,'-o')
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2','y_3','Location','NorthWest');

Figure contains an axes object. The axes object with xlabel Time t, ylabel Solution y contains 3 objects of type line. These objects represent y_1, y_2, y_3.

ローカル関数

ここでは、DDE ソルバー dde23 が解を計算するために呼び出すローカル補助関数を紹介しています。あるいは、これらの関数を独自のファイルとして MATLAB パスのディレクトリに保存することもできます。

function dydt = ddefun(t,y,Z) % equation being solved
  ylag1 = Z(:,1);
  ylag2 = Z(:,2);

  dydt = [ylag1(1); 
          ylag1(1)+ylag2(2); 
          y(2)];
end
%-------------------------------------------
function s = history(t) % history function for t <= 0
  s = ones(3,1);
end
%-------------------------------------------

参考

dde23 | ddesd | ddensd | deval