# ddensd
求解中立型时滞微分方程 (DDE)
函数库: TyMath
# 语法
sol = ddensd(ddefun,dely,delyp,history,tspan)
sol = ddensd(ddefun,dely,delyp,history,tspan,options)
# 说明
sol = ddensd(ddefun ,dely,delyp,history,tspan) 计算以下形式的中立型时滞微分方程组的积分：
，其中
t 是表示时间的自变量；
dyi 是任何 p 解时滞；
dypi 是任何 q 导数时滞。
sol = ddensd(ddefun ,dely,delyp,history,tspan,options) 将默认积分参数替换为在使用 ddeset 函数创建的结构体 options 中指定的参数。
# 示例
带有两个时滞的中立型 DDE
对 Paul 提出的以下中立型 DDE 求解
在编辑器中打开一个新程序文件。此文件包含一个主函数和四个局部函数。
将一阶 DDE 定义为名为 ddefun 的局部函数。
using TyMath
using TyPlot
function ddefun(t,y,ydel,ypdel)
yp = @. 1 + y - 2*ydel^2 - ypdel;
end
将解时滞定义为名为 dely 的局部函数。
function dely(t,y)
dy = t/2;
end
将导数时滞定义为名为 delyp 的局部函数。
function delyp(t,y)
dyp = t-pi;
end
将解历史记录定义为名为 history 的局部函数。
function history(t)
y = cos(t);
end
定义积分区间并使用 ddensd 对 DDE 求解。将此代码添加到主函数。
tspan = [0,pi];
sol = ddensd(ddefun,dely,delyp,history,tspan);
在
tn = LinRange(0,pi,100);
yn = interp1(sol.x,sol.y',tn);
绘制结果。将此代码添加到主函数。
plot(tn,yn);
xlim([0,pi]);
ylim([-1.2,1.2]);
xlabel("time t");
ylabel("solution y");
# 输入参数
ddefun - 导函数函数句柄
导函数，指定为函数句柄，其语法为 yp = ddefun(t,y,ydel,ypdel)。ddefun 的参数见下表所述。
| ddefun 参数 | 说明 |
|---|---|
| t | 表示当前时间值 t 的标量值。 |
| y | 以公式 1 表示 y(t) 的向量。此向量的大小为 n x 1，其中 n 是您需要求解的方程组中的方程数。 |
| ydel | 由表示 y[dyi] 的列 ydel[:,i] 构成的矩阵。此矩阵的大小为 n x p，其中 n 是您需要求解的方程组中的方程数，p 是公式 1 中的 y[dy] 项的数量。 |
| ypdel | 由表示 y[dypj] 的列 ypdel[:,j] 构成的矩阵。此矩阵的大小为 n x q，其中 n 是您需要求解的方程组中的方程数，q 是 公式 1 中的 y '[dyp] 项的数量。 |
| yp | ddefun 返回的结果。它是一个 n x 1 向量，其元素表示公式 1 的右侧。 |
dely - 解时滞函数句柄 | 向量
解时滞，指定为函数句柄，会在公式 1 中返回 dy1,..., dyp。您也可以通过向量形式传递固定时滞。
如果将 dely 指定为函数句柄，语法必须为 dy = dely(t,y)。此函数的参数见下表所述。
| dely 参数 | 说明 |
|---|---|
| t | 表示当前时间值 t 的标量值。 |
| y | 以公式 1 表示 y(t) 的向量。此向量的大小为 n x 1，其中 n 是您需要求解的方程组中的方程数。 |
| dy | dely 函数返回的向量，其值为公式 1 中的解时滞 dyi。此向量的大小为 p x 1，其中 p 是方程中的解时滞数。每个元素必须小于或等于 t。 |
如果想要指定 dyi = t – τi 形式的固定解时滞，则 dely 必须是向量，其中 dely[i] = τi。此向量中的每个值必须大于或等于零。
如果问题中不存在 dy，则将 dely 设置为 nothing。
数据类型： Function | single | double
delyp - 导数时滞函数句柄 | 向量
导数时滞，指定为函数句柄，会在公式 1 中返回 dyp1,..., dypq。您也可以通过向量形式传递固定时滞。
如果 delyp 是一个函数句柄，则其语法必须是 dyp = delyp(t,y)。此函数的参数见下表所述。
| dely 参数 | 说明 |
|---|---|
| t | 表示当前时间值 t 的标量值。 |
| y | 以公式 1 表示 y(t) 的向量。此向量的大小为 n x 1，其中 n 是您需要求解的方程组中的方程数。 |
| dyp | delyp 函数返回的向量，其值为公式 1 中的导数时滞 dypj。此向量的大小必须为 q x 1，其中 q 是方程中的解时滞 dypj 的数量。dyp 的每个元素必须小于 t。此限制有一个例外：如果您是在求解初始值 DDE，则 dyp 的值在 t = t0 时可以等于 t。 |
如果想要指定 dypi = t – τi 形式的固定导数时滞，则 delyp 必须是向量，其中 delyp(j) = τj。此向量中的每个值必须大于零。对于这一限制，当求解中立型 DDE 初始值时，会出现例外情况。在这种情况下，当 t = t0 时，delyp 可以等于零。
如果问题中不存在 dyp，则将 delyp 设置为 nothing。
数据类型： Function | AbstractFloat
history - 解历史记录函数句柄 | 列向量 | 结构体（sol，来自之前的积分） | 1 x 2 元组
解历史记录，指定为函数句柄、列向量、sol 结构体（来自之前积分）或元胞数组。这是 t ≤ t0 处的解。
如果历史记录随着时间发生变化，则将解历史记录指定为一个语法为 y = history(t) 的函数句柄。此函数为 t <= t0 返回由 y(t) 的近似解组成的 n x 1 向量。此向量的长度为 n，是您需要求解的方程组中的方程数；
如果 y(t) 是常量，您可以将 history 指定为一个 n-1 的常量值向量；
如果想要调用 ddensd 以将之前的积分继续到 t0，可以将历史记录指定为来自之前积分的输出 sol；
如果要求解初始值 DDE，请将历史记录指定为元胞数组 {y0, yp0}。第一个元素 y0 是初始值 y(t0) 的列向量。第二个元素 yp0 是元素为初始导数 y '(t0) 的列向量。这些向量必须保持一致，也就是说它们必须满足 t0 处的公式 1。
数据类型： Function | AbstractFloat | struct | Tuple
tspan - 积分区间1 x 2 向量
积分区间，指定为向量 [t0 tf]。第一个元素 t0 是 t 的初始值。第二个元素 tf 是 t 的最终值。t0 的值必须小于 tf。
数据类型： Number
options - 可选积分参数ddeset 返回的结构体
可选积分参数，指定为由 ddeset 函数创建和返回的结构体。一些常用属性包括：'RelTol'、'AbsTol' 和 'Events'。
# 输出参数
sol - 解结构体
解，以包含下列字段的结构体形式返回。
| 参数 | 说明 |
|---|---|
| sol.x | ddensd 选择的网格。 |
| sol.y | 网格点处的 y(t) 近似值。 |
| sol.yp | 网格点处的 y′(t) 近似值。 |
| sol.solver | 标识求解器 'ddensd' 的字符向量。 |
# 详细信息
初始值中立型时滞微分方程
初始值 DDE 具有 dyi≥t0 并且 dypj≥t0（对于所有 i 和 j）。在 t = t0 处，所有时滞项降低为 y[dyi] = y[t0] 和 y '[dypj] = y '[t0]：
对于 t > t0，所有导数时滞必须满足 dyp < t。
求解初始值中立型 DDE 时，必须为 ddensd 提供 y '(t0)。为此，请将 history 指定为元胞数组 (Y0,YP0)。其中，Y0 是初始值 y(t0) 的列向量，YP0 是初始导数 y '(t0) 的列向量。
# 参考文献
[1] Paul, C.A.H. “A Test Set of Functional Differential Equations.” Numerical Analysis Reports. No. 243. Manchester, UK: Math Department, University of Manchester, 1994.
[2] Shampine, L.F. “Dissipative Approximations to Neutral DDEs.” Applied Mathematics & Computation. Vol. 203, Number 2, 2008, pp. 641–648.