MATLAB 数学应用 微分方程 时滞微分方程 中立型的初始值DDE

本文讲述了如何使用 ddensd 求解具有时间相关时滞的初始值 DDE（时滞微分方程）方程组。

方程是：

$y^{\prime }(t)=2cos(2t)y(\frac{t}{2}{)}^{2cos(t)}+log(y^{\prime }(\frac{t}{2}))-log(2cos(t))-sin(t).$

此方程是初始值 DDE，因为在 $t_0$ 处时滞为零。因此，不需要历史解来计算解，只需要初始值：

$y(0)=1,$

$y^{\prime }(0)=s.$

$s$ 是 $2+log(s)-log(2)=0$ 的解。满足此方程的 $s$ 的值是 $s_1=2$ 和 $s_2=0.4063757399599599$。

由于方程中的时滞存在于 $y^{\prime }$ 项中，因此该方程称为中立型 DDE。

要在 MATLAB 中求解此方程，您需要先编写方程和时滞的代码，然后调用时滞微分方程求解器 ddensd，后者是中立型方程的求解器。您可以将所需的函数作为局部函数包含在文件末尾，或者将它们作为单独的文件保存在 MATLAB 路径上的目录中。

编写时滞代码

首先，编写一个匿名函数来定义方程中的迟滞。由于 y 和 y′ 都 $t_2$ 形式的迟滞，因此只需要一个函数定义。此时滞函数后来传递给求解器两次，一次表示 $y$ 的时滞，一次表示 $y^{\prime }$ 的时滞。

delay = @(t,y) t/2; 

编写方程代码

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

• t 是时间（自变量）。
• y 是解（因变量）。
• ydel 包含 y 的时滞。
• ypdel 包含 y 的时滞。

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

• ydel→$y(\frac{t}{2})$
• ypdel →$y^{\prime }(\frac{t}{2})$

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

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

求解方程

最后，定义积分区间 [$t_0$ $t_f$] 和初始值，然后使用 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);

对解进行绘图

解结构体 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');

局部函数

此处列出了 DDE 求解器 ddensd 为计算解而调用的局部辅助函数。您也可以将这些函数作为它们自己的文件保存在 MATLAB 路径上的目录中。

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