大家好，又见面了，我是你们的朋友全栈君。如果您正在找激活码,请点击查看最新教程,关注关注公众号 “全栈程序员社区” 获取激活教程,可能之前旧版本教程已经失效.最新Idea2022.1教程亲测有效,一键激活。
Jetbrains全系列IDE稳定放心使用
本文讲述了如何使用 dde23 对具有常时滞的DDE（时滞微分方程）方程组求解。
方程组为：
y 1 ′ ( t ) = y 1 ( t − 1 ) y’_1(t)=y_1(t−1) y1′(t)=y1(t−1)
y 2 ′ ( t ) = y 1 ( t − 1 ) + y 2 ( t − 0.2 ) y’_2(t)=y_1(t-1)+y_2(t-0.2) y2′(t)=y1(t−1)+y2(t−0.2)
y 3 ′ ( t ) = y 2 ( t ) y’_3(t)=y_2(t) y3′(t)=y2(t).
t≤0 的历史解函数是常量 y 1 ( t ) = y 2 ( t ) = y 3 ( t ) = 1 y_1(t)=y_2(t)=y_3(t)=1 y1(t)=y2(t)=y3(t)=1。
方程中的时滞仅存在于 y 项中，并且时滞本身是常量，因此各方程构成常时滞方程组。
要在 MATLAB 中求解此方程组，您需要先编写方程组、时滞和历史解的代码，然后再调用时滞微分方程求解器 dde23，该求解器适用于具有常时滞的方程组。您可以将所需的函数作为局部函数包含在文件末尾，或者将它们作为单独的命名文件保存在 MATLAB 路径上的目录中。
编写时滞代码
首先，创建一个向量来定义方程组中的时滞。此方程组有两种不同时滞：
- 在第一个分量 y 1 ( t − 1 ) y_1(t−1) y1(t−1) 中时滞为 1。
- 在第二个分量 y 2 ( t − 0.2 ) y_2(t−0.2) y2(t−0.2) 中时滞为 0.2。
dde23 接受时滞的向量参数，其中每个元素是一个分量的常时滞。
lags = [1 0.2];
编写方程代码
现在，创建一个函数来编写方程的代码。此函数应具有签名 dydt = ddefun(t,y,Z)，其中：
t 是时间（自变量）。
y 是解（因变量）。
Z(:,j) 用于逼近时滞 y ( t − τ j ) y(t−τ_j) y(t−τj)，其中常时滞 τ j τ_j τj 由 lags(j) 给定。
求解器会自动将这些输入传递给该函数，但是变量名称决定如何编写方程代码。在这种情况下：
- Z ( : , 1 ) → y 1 ( t − 1 ) Z(:,1) → y_1(t−1) Z(:,1) → y1(t−1)
- Z ( : , 2 ) → y 2 ( t − 0.2 ) Z(:,2) → y_2(t−0.2) 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 ≤ t 0 t≤t_0 t≤t0 的解。
function s = history(t)
s = ones(3,1);
end
求解方程
最后，定义积分区间 [ t 0 t f ] [t_0 t_f] [t0 tf] 并使用 dde23 求解器对 DDE 求解。
tspan = [0 5];
sol = dde23(@ddefun, lags, @history, tspan);
对解进行绘图
解结构体 sol 具有字段 sol.x 和 sol.y，这两个字段包含求解器在这些时间点所用的内部时间步和对应的解。（如果您需要在特定点的解，可以使用 deval 来计算在特定点的解。）
绘制三个解分量对时间的图。
plot(sol.x,sol.y,'-o')
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2','y_3','Location','NorthWest');
局部函数
此处列出了 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
%-------------------------------------------
发布者：全栈程序员-站长，转载请注明出处：https://javaforall.net/188171.html原文链接：https://javaforall.net
