matlab使用教程(34)—求解时滞微分方程(2)

1.具有状态依赖时滞的 DDE

        以下示例说明如何使用 ddesd 对具有状态依赖时滞的 DDE（时滞微分方程）方程组求解。Enright 和Hayashi [1] 将此 DDE 方程组用作测试问题。方程组为：

        方程中的时滞仅出现在 y 项中。时滞仅取决于第二个分量 y 2 t 的状态，因此这些方程构成状态依赖时滞方程组。

        要在 MATLAB® 中求解此方程组，您需要先编写方程组、时滞和历史解的代码，然后再调用时滞微分方程求解器 ddesd，该求解器适用于具有状态依赖时滞的方程组。您可以将所需的函数作为局部函数包含在文件末尾（如本处所示），或者将它们作为单独的命名文件保存在 MATLAB 路径上的目录中。

1.1编写时滞代码

function d = dely(t,y)
d = exp(1 - y(2));
end

1.2编写方程代码

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

function dydt = ddefun(t,y,Z)
dydt = [y(2);
 -Z(2,1)*y(2)^2*exp(1 - y(2))];
end

1.3编写历史解代码

        接下来，创建一个函数来定义历史解。历史解是时间 t ≤ t 0 的解。

function v = history(t) % history function for t < t0
v = [log(t); 
 1./t];
end

1.4求解方程

        最后，定义积分区间并使用 ddesd 求解器对 DDE 求解。

tspan = [0.1 5];
sol = ddesd(@ddefun, @dely, @history, tspan);

1.5对解进行绘图

        解结构体 sol 具有字段 sol.x 和 sol.y，这两个字段包含求解器在这些时间点所用的内部时间步和对应的解。（如果您需要在特定点的解，可以使用 deval 来计算在特定点的解。）使用历史解函数绘制两个解分量对时间的图，以计算积分区间内的解析解来进行比较。

ta = linspace(0.1,5);
ya = history(ta);
plot(ta,ya,sol.x,sol.y,'o')
legend('y_1 exact','y_2 exact','y_1 ddesd','y_2 ddesd')
xlabel('Time t')
ylabel('Solution y')
title('D1 Problem of Enright and Hayashi')

1.6局部函数

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

function dydt = ddefun(t,y,Z) % equation being solved
dydt = [y(2); 
 -Z(2,1).*y(2)^2.*exp(1 - y(2))];
end
%-------------------------------------------
function d = dely(t,y) % delay for y
d = exp(1 - y(2));
end
%-------------------------------------------
function v = history(t) % history function for t < t0
v = [log(t); 
 1./t];
end

2具有不连续性的心血管模型 DDE

        此示例说明如何使用 dde23 对具有不连续导数的心血管模型求解。此示例最初由 Ottesen [1] 提出。方程组为：

        该方程组受外周压的巨大影响，外周压会从 R = 1 . 05 急剧减少到 R = 0 . 84 ，从 t = 600 处开始。因此，该方程组在 t = 600 处的低阶导数具有不连续性。常历史解由以下物理参数定义

        要在 MATLAB® 中求解此方程组，您需要先编写方程组、参数、时滞和历史解的代码，然后再调用时滞微分方程求解器 dde23，该求解器适用于具有常时滞的方程组。您可以将所需的函数作为局部函数包含在文件末尾（如本处所示），或者将它们作为单独的命名文件保存在 MATLAB 路径上的目录中。

2.1定义物理参数

首先，将问题的物理参数定义为结构体中的字段。

p.ca = 1.55;
p.cv = 519;
p.R = 1.05;
p.r = 0.068;
p.Vstr = 67.9;
p.alpha0 = 93;
p.alphas = 93;
p.alphap = 93;
p.alphaH = 0.84;
p.beta0 = 7;
p.betas = 7;
p.betap = 7;
p.betaH = 1.17;
p.gammaH = 0;

tau = 4;

2.2编写方程代码

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

        求解器自动将前三个输入传递给函数，变量名称决定如何编写方程代码。调用求解器时，参数结构体 p 将传递给函数。在本例中，时滞表示为：

function dydt = ddefun(t,y,Z,p)
 if t <= 600
 p.R = 1.05;
 else
 p.R = 0.21 * exp(600-t) + 0.84;
 end 
 ylag = Z(:,1);
 Patau = ylag(1);
 Paoft = y(1);
 Pvoft = y(2);
 Hoft = y(3);
 dPadt = - (1 / (p.ca * p.R)) * Paoft ...
 + (1/(p.ca * p.R)) * Pvoft ...
 + (1/p.ca) * p.Vstr * Hoft;
 dPvdt = (1 / (p.cv * p.R)) * Paoft...
 - ( 1 / (p.cv * p.R)...
 + 1 / (p.cv * p.r) ) * Pvoft;
 Ts = 1 / ( 1 + (Patau / p.alphas)^p.betas );
 Tp = 1 / ( 1 + (p.alphap / Paoft)^p.betap );
 dHdt = (p.alphaH * Ts) / (1 + p.gammaH * Tp) ...
 - p.betaH * Tp;
 dydt = [dPadt; dPvdt; dHdt];
end

2.3编写历史解代码

接下来，创建一个向量来定义三个分量 P a 、 P v 和 H 的常历史解。历史解是时间 t ≤ t 0 的解。

P0 = 93;
Paval = P0;
Pvval = (1 / (1 + p.R/p.r)) * P0;
Hval = (1 / (p.R * p.Vstr)) * (1 / (1 + p.r/p.R)) * P0;
history = [Paval; Pvval; Hval];

2.4 求解方程

        使用 ddeset 来指定在 t = 600 处存在不连续性。最后，定义积分区间 并使用 dde23 求解器对 DDE 求解。使用匿名函数指定 ddefun 以传入参数结构体 p 。

options = ddeset('Jumps',600);
tspan = [0 1000];
sol = dde23(@(t,y,Z) ddefun(t,y,Z,p), tau, history, tspan, options);

2.4对解进行绘图

        解结构体 sol 具有字段 sol.x 和 sol.y，这两个字段包含求解器在这些时间点所用的内部时间步和对应的解。（如果您需要在特定点的解，可以使用 deval 来计算在特定点的解。）绘制第三个解分量（心率）对时间的图。

plot(sol.x,sol.y(3,:))
title('Heart Rate for Baroreflex-Feedback Mechanism.')
xlabel('Time t')
ylabel('H(t)')

2.5局部函数

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

function dydt = ddefun(t,y,Z,p) % equation being solved
 if t <= 600
 p.R = 1.05;
 else
 p.R = 0.21 * exp(600-t) + 0.84;
 end 
 ylag = Z(:,1);
 Patau = ylag(1);
 Paoft = y(1);
 Pvoft = y(2);
 Hoft = y(3);
 dPadt = - (1 / (p.ca * p.R)) * Paoft ...
 + (1/(p.ca * p.R)) * Pvoft ...
 + (1/p.ca) * p.Vstr * Hoft;
 dPvdt = (1 / (p.cv * p.R)) * Paoft...
 - ( 1 / (p.cv * p.R)...
 + 1 / (p.cv * p.r) ) * Pvoft;
 Ts = 1 / ( 1 + (Patau / p.alphas)^p.betas );
 Tp = 1 / ( 1 + (p.alphap / Paoft)^p.betap );
 dHdt = (p.alphaH * Ts) / (1 + p.gammaH * Tp) ...
 - p.betaH * Tp;
 dydt = [dPadt; dPvdt; dHdt];
end