1中立型 DDE
以下示例说明如何使用 ddensd 求解中立型 DDE(时滞微分方程),其中时滞出现在导数项中。此问题最初由 Paul [1] 提出。方程是:
由于该方程在 y
′ 项中存在时滞,因此该方程称为中立型 DDE。如果时滞仅出现在 y 项中,则根据时滞的形式,方程将是常时滞或状态依赖 DDE。要在 MATLAB® 中求解此方程,您需要先编写方程、时滞和历史解的代码,然后再调用时滞微分方程求解器
ddensd。您可以将这些作为局部函数包含在文件末尾(如本处所示),或者将它们作为单独的文件保存在 MATLAB 路径上的目录中。
1.1编写时滞代码
function dy = dely(t,y)
dy = t/2;
end
function dyp = delyp(t,y)
dyp = t-pi;
end
在此示例中,y
和
y
′ 分别仅有一个时滞。如果有更多时滞,则您可以将它们添加到这些相同的函数文件中,这样函数将返回向量而不是标量。
1.2编写方程代码
现在,创建一个函数来编写方程代码。此函数应具有签名
yp = ddefun(t,y,ydel,ypdel)
,其中:
function yp = ddefun(t,y,ydel,ypdel)
yp = 1 + y - 2*ydel^2 - ypdel;
end
1.3编写历史解代码
接下来,创建一个函数来定义历史解。历史解是时间 t
≤
t
0
的解。
function y = history(t)
y = cos(t);
end
1.4 求解方程
最后,定义积分区间
并使用
ddensd
求解器对 DDE 求解。
tspan = [0 pi];
sol = ddensd(@ddefun, @dely, @delyp, @history, [0,pi]);
1.5对解进行绘图
解结构体 sol
具有字段
sol.x
和
sol.y,这两个字段包含求解器在这些时间点所用的内部时间步和对应的解。但是,您可以使用
deval 计算在特定点的解。在
0
和
pi 之间以 20 个等间距点计算解。绘制计算解和历史解对解析解的图。
tn = linspace(0,pi,20);
yn = deval(sol,tn);
th = linspace(-pi,0);
yh = history(th);
ta = linspace(0,pi);
ya = cos(ta);
plot(th,yh,tn,yn,'o',ta,ya)
legend('History','Numerical','Analytical','Location','NorthWest')
xlabel('Time t')
ylabel('Solution y')
title('Example of Paul with 1 Equation and 2 Delay Functions')
axis([-3.5 3.5 -1.5 1.5])
1.6局部函数
此处列出了 DDE 求解器 ddensd 为计算解而调用的局部辅助函数。您也可以将这些函数作为它们自己的文件保存在 MATLAB 路径上的目录中。
function yp = ddefun(t,y,ydel,ypdel) % equation being solved
yp = 1 + y - 2*ydel^2 - ypdel;
end
%-------------------------------------------
function dy = dely(t,y) % delay for y
dy = t/2;
end
%-------------------------------------------
function dyp = delyp(t,y) % delay for y'
dyp = t-pi;
end
%-------------------------------------------
function y = history(t) % history function for t < 0
y = cos(t);
end
2中立型的初始值 DDE
以下示例说明如何使用 ddensd 求解具有时间依赖时滞的初始值 DDE(时滞微分方程)方程组。此示例最初由 Jackiewicz [1] 提出。
由于方程中的时滞存在于 y
′
项中,因此该方程称为中立型 DDE。
要在 MATLAB® 中求解此方程,您需要先编写方程和时滞的代码,然后调用时滞微分方程求解器ddensd,后者是中立型方程的求解器。您可以将所需的函数作为局部函数包含在文件末尾(如本处所示),或者将它们作为单独的文件保存在 MATLAB 路径上的目录中。
2.1编写时滞代码
首先,编写一个匿名函数来定义方程中的迟滞。由于 y
和
y′ 都有 t/2 形式的迟滞,因此只需要一个函数定 义。此时滞函数后来传递给求解器两次,一次表示
y
的时滞,一次表示
y
′
的时滞。
delay = @(t,y) t/2;
2.2编写方程代码
现在,创建一个函数来编写方程代码。此函数应具有签名 yp = ddefun(t,y,ydel,ypdel)
,其中:
function yp = ddefun(t,y,ydel,ypdel)
yp = 2*cos(2*t)*ydel^(2*cos(t)) + log(ypdel) - log(2*cos(t)) - sin(t);
end
2.3求解方程
最后,定义积分区间
和初始值,然后使用
ddensd 求解器求解 DDE。通过在第四个输入参数的元胞数组中指定初始值,将初始值传递给求解器。
tspan = [0 0.1];
y0 = 1;
s1 = 2;
sol1 = ddensd(@ddefun, delay, delay, {y0,s1}, tspan);
第二次求解方程,这次使用
s
的备选值作为初始条件。
s2 = 0.4063757399599599;
sol2 = ddensd(@ddefun, delay, delay, {y0,s2}, tspan);
2.4对解进行绘图
解结构体 sol1
和
sol2
具有字段
x
和
y,这些字段包含求解器在这些时间点所用的内部时间步和对应的解。但是,您可以使用
deval 计算在特定点的解。绘制两个解以比较结果。
plot(sol1.x,sol1.y,sol2.x,sol2.y);
legend('y''(0) = 2','y''(0) = .40637..','Location','NorthWest');
xlabel('Time t');
ylabel('Solution y');
title('Two Solutions of Jackiewicz''s Initial-Value NDDE');