matlab使用教程(35)—求解时滞微分方程(3)

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