具有状态依赖时滞的 DDE

打开实时脚本

以下示例说明如何使用 ddesd 对具有状态依赖时滞的 DDE（时滞微分方程）方程组求解。Enright 和 Hayashi [1] 将此 DDE 方程组用作测试问题。

方程组为：

${y_{1}}^{'} (t) = y_{2} (t),$

${y_{2}}^{'} (t) = - y_{2} (e^{1 - y_{2} (t)}) \cdot y_{2} {(t)}^{2} \cdot e^{1 - y_{2} (t)} .$

$t \leq 0.1$ 的历史解函数是解析解

$y_{1} (t) = \log (t),$

$y_{2} (t) = \frac{1}{t} .$

方程中的时滞仅出现在 $y$ 项中。时滞仅取决于第二个分量 $y_{2} (t)$ 的状态，因此这些方程构成状态依赖时滞方程组。

要在 MATLAB® 中求解此方程组，您需要先编写方程组、时滞和历史解的代码，然后再调用时滞微分方程求解器 ddesd，该求解器适用于具有状态依赖时滞的方程组。您可以将所需的函数作为局部函数包含在文件末尾（如本处所示），或者将它们作为单独的命名文件保存在 MATLAB 路径上的目录中。

编写时滞代码

首先，编写一个函数来定义方程组中的时滞。此方程组中唯一存在的时滞位于项 $- y_{2} (e^{1 - y_{2} (t)})$ 中。

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

注意：所有函数都作为局部函数包含在示例的末尾。

编写方程代码

现在，创建一个函数来编写方程的代码。此函数应具有签名 dydt = ddefun(t,y,Z)，其中：

t 是时间（自变量）。
y 是解（因变量）。
Z(n,j) 对时滞 $y_{n} (d (j))$ 求近似值，其中时滞 $d (j)$ 由 dely(t,y) 的分量 j 给出。

求解器会自动将这些输入传递给该函数，但是变量名称决定如何编写方程代码。在这种情况下：

Z(2,1) $\to y_{2} (e^{1 - y_{2} (t)})$

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

编写历史解代码

接下来，创建一个函数来定义历史解。历史解是时间 $t \leq t_{0}$ 的解。

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

求解方程

最后，定义积分区间 $[t_{0} t_{f}]$ 并使用 ddesd 求解器对 DDE 求解。

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

对解进行绘图

解结构体 sol 具有字段 sol.x 和 sol.y，这两个字段包含求解器在这些时间点所用的内部时间步和对应的解。（如果您需要在特定点的解，可以使用 deval 来计算在特定点的解。）

使用历史解函数绘制两个解分量对时间的图，以计算积分区间内的解析解来进行比较。

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')

Figure contains an axes object. The axes object with title D1 Problem of Enright and Hayashi, xlabel Time t, ylabel Solution y contains 4 objects of type line. One or more of the lines displays its values using only markers These objects represent y_1 exact, y_2 exact, y_1 ddesd, y_2 ddesd.

局部函数

此处列出了 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