状態依存の遅延をもつ DDE

この例では、ddesd を使用して、状態依存の遅延を伴う DDE (遅延微分方程式) 系の解を求める方法を説明します。この DDE 系は、Enright と Hayashi によるテスト問題として使用されました [1]。

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

$${{\mathit{y}}_1}^{\prime }\left(\mathit{t}\right)={\mathit{y}}_2\left(\mathit{t}\right),$$

$${{\mathit{y}}_2}^{\prime }\left(\mathit{t}\right)=-{\mathit{y}}_2\left({\mathit{e}}^{1-{\mathit{y}}_2\left(\mathit{t}\right)}\right)\cdot {\mathit{y}}_2{\left(\mathit{t}\right)}^2\cdot {\mathit{e}}^{1-{\mathit{y}}_2\left(\mathit{t}\right)}.$$

$\mathit{t}\le 0.1$ での履歴関数は次の解析解です。

$${\mathit{y}}_1\left(\mathit{t}\right)=\log \left(\mathit{t}\right),$$

$${\mathit{y}}_2\left(\mathit{t}\right)=\frac{1}{\mathit{t}}.$$

方程式の時間遅延は $\mathit{y}$ の項にのみ存在します。遅延は 2 番目の要素 ${\mathit{y}}_2\left(\mathit{t}\right)$ の状態にのみ依存するため、方程式は "状態依存の遅延" の方程式系を形成します。

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

遅延のコード化

最初に、方程式系の時間遅延を定義する関数を作成します。この方程式系に存在する唯一の遅延は項 $-{\mathit{y}}_2\left({\mathit{e}}^{1-{\mathit{y}}_2\left(\mathit{t}\right)}\right)$ にあります。

function d = dely(t,y)
d = exp(1 - y(2));
end

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

方程式のコード化

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

• t は時間 (独立変数) です。

• y は解 (従属変数) です。

• Z(n,j) は遅延 ${\mathit{y}}_{\mathit{n}}\left(\mathit{d}\left(\mathit{j}\right)\right)$ を近似します。ここで、遅延 $\mathit{d}\left(\mathit{j}\right)$ は dely(t,y) の要素 j によって与えられます。

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

• Z(2,1)$\to {\mathit{y}}_2\left({\mathit{e}}^{1-{\mathit{y}}_2\left(\mathit{t}\right)}\right)$

function dydt = ddefun(t,y,Z)
dydt = [y(2);
       -Z(2,1)*y(2)^2*exp(1 - y(2))];
end

解の履歴のコード化

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

function v = history(t) % history function for t < t0
v = [log(t); 
     1./t];
end

方程式の求解

最後に、積分区間 $\left[{\mathit{t}}_0{\mathit{t}}_{\mathit{f}}\right]$ を定義し、ddesd ソルバーを使用して DDE を解きます。

tspan = [0.1 5];
sol = ddesd(@ddefun, @dely, @history, tspan);

解のプロット

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

履歴関数を使用して 2 つの解要素を時間に対してプロットし、比較のために積分区間内の解析解を計算します。

ta = linspace(0.1,5);
ya = history(ta);

plot(ta,ya,sol.x,sol.y,'o')
legend('y_1 exact','y_2 exact','y_1 ddesd','y_2 ddesd')
xlabel('Time t')
ylabel('Solution y')
title('D1 Problem of Enright and Hayashi')

ローカル関数

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

function dydt = ddefun(t,y,Z) % equation being solved
dydt = [y(2); 
       -Z(2,1).*y(2)^2.*exp(1 - y(2))];
end
%-------------------------------------------
function d = dely(t,y) % delay for y
d = exp(1 - y(2));
end
%-------------------------------------------
function v = history(t) % history function for t < t0
v = [log(t); 
     1./t];
end
%-------------------------------------------

参照

[1] Enright, W.H. and H. Hayashi."The Evaluation of Numerical Software for Delay Differential Equations."In Proceedings of the IFIP TC2/WG2.5 working conference on Quality of numerical software: assessment and enhancement.(R.F. Boisvert, ed.).London, UK:Chapman & Hall, Ltd., pp. 179-193.

参考

ddensd | ddesd | dde23 | deval

関連するトピック

• 遅延微分方程式の求解
• 不連続をもつ心臓血管モデル DDE