matlab用dde23求解带有固定时滞的时滞微分方程

一个同学咨询的带有固定时滞的时滞微分方程求解，故分享一下matlab中dde23的用法

dde23函数调用方法

sol = dde23(ddefun,lags,history,tspan,options)

dde23 跟踪不连续性并使用显式 Runge-Kutta (2,3) 对和插值对 ode23 求积分。它通过迭代来采用超过时滞的步长。

举例：

t≤0 的历史解函数是常量 y1(t)=y2(t)=y3(t)=1。

方程中的时滞仅存在于 y 项中，并且时滞本身是常量，因此各方程构成常时滞方程组。

要在 MATLAB 中求解此方程组，需要先编写方程组、时滞和历史解的代码，然后再调用时滞微分方程求解器 dde23，该求解器适用于具有常时滞的方程组。可以将所需的函数作为局部函数或者将它们作为单独的命名文件保存在 MATLAB 路径上的目录中。

编写时滞代码

首先，创建一个向量来定义方程组中的时滞。此方程组有两种不同时滞：

• 在第一个分量 y1(t−1) 中时滞为 1。
• 在第二个分量 y2(t−0.2) 中时滞为 0.2。

dde23 接受时滞的向量参数，其中每个元素是一个分量的常时滞。

lags = [1 0.2];

编写方程代码

现在，创建一个函数来编写方程的代码。此函数应变换为这种格式：

dydt = ddefun(t,y,Z)

其中：

• t 是时间（自变量）。
• y 是解（因变量）。
• Z(:,j) 用于逼近时滞 y(t−τj)，其中常时滞 τj 由 lags(j) 给定。

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

• Z(:,1) → y1(t−1)
• Z(:,2) → y2(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≤t0 的解。

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

求解方程

最后，定义积分区间 [t0 tf] 并使用 dde23 求解器对 DDE 求解。

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

对解进行绘图

求解的结构体 sol 具有字段 sol.x 和 sol.y，这两个字段包含求解器在这些时间点所用的内部时间步和对应的解。绘制三个解分量对时间的图。

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