具有不连续性的心血管模型 DDE
此示例说明如何使用 dde23
对具有不连续导数的心血管模型求解。此示例最初由 Ottesen [1] 提出。
方程组为:
和 的项分别是同一方程在有时滞和没有时滞状态下的变体。 和 分别代表在有时滞和没有时滞状态下的平均动脉压。
此问题有许多物理参数:
动脉顺应性
静脉顺应性
外周阻力
静脉流出阻力
心博量
典型平均动脉压
该方程组受外周压的巨大影响,外周压会从 急剧减少到 ,从 处开始。因此,该方程组在 处的低阶导数具有不连续性。
常历史解由以下物理参数定义
要在 MATLAB® 中求解此方程组,您需要先编写方程组、参数、时滞和历史解的代码,然后再调用时滞微分方程求解器 dde23
,该求解器适用于具有常时滞的方程组。您可以将所需的函数作为局部函数包含在文件末尾(如本处所示),或者将它们作为单独的命名文件保存在 MATLAB 路径上的目录中。
定义物理参数
首先,将问题的物理参数定义为结构体中的字段。
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
来表示项 的方程中的常时滞 。
tau = 4;
编写方程代码
现在,创建一个函数来编写方程的代码。此函数应具有签名 dydt = ddefun(t,y,Z,p)
,其中:
t
是时间(自变量)。y
是解(因变量)。Z(n,j)
对时滞 求近似值,其中时滞 由dely(t,y)
的分量j
给出。p
是可选的第四个输入,用于传入参数值。
求解器自动将前三个输入传递给函数,变量名称决定如何编写方程代码。调用求解器时,参数结构体 p
将传递给函数。在本例中,时滞表示为:
Z(:,1)
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
注意:所有函数都作为局部函数包含在示例的末尾。
编写历史解代码
接下来,创建一个向量来定义三个分量 、 和 的常历史解。历史解是时间 的解。
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];
求解方程
使用 ddeset
来指定在 处存在不连续性。最后,定义积分区间 并使用 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);
对解进行绘图
解结构体 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)')
局部函数
此处列出了 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
参考
[1] Ottesen, J. T.“Modelling of the Baroreflex-Feedback Mechanism with Time-Delay.”J. Math.Biol.Vol. 36, Number 1, 1997, pp. 41–63.
另请参阅
ddensd
| ddesd
| dde23
| deval