不连续的心血管模型DDE

这个例子展示了如何使用dde23解一个具有不连续导数的心血管模型。该示例最初由Ottesen[1]提出。

方程组是

$\dot{P_{一个}} （ t ）＝ - \frac{1}{c_{一个} R} P_{一个} （ t ） + \frac{1}{c_{一个} R} P_{v} （ t ） + \frac{1}{c_{一个}} V_{str} （ P_{一个}^{τ} （ t ）） H （ t ）$

${\dot{P}}_{v} （ t ）＝ \frac{1}{c_{v} R} P_{一个} （ t ） - （ \frac{1}{c_{v} R} + \frac{1}{c_{v} r} ） P_{v} （ t ）$

$\dot{H} （ t ）＝ \frac{α_{H} T_{年代}}{1 + γ_{H} T_{p}} - β_{H} T_{p} ．$

的条款 $T_{年代}$ 而且 $T_{p}$ 是同一方程的有或没有时滞的变化形式。 $P_{一个}^{τ}$ 而且 $P_{一个}$ 分别表示有延时和无延时的平均动脉压。

$T_{年代} ＝ \frac{1}{1 + {（ \frac{P_{一个}^{τ}}{α_{年代}} ）}^{β_{年代}}}$

$T_{p} ＝ \frac{1}{1 + {（ \frac{P_{一个}}{α_{p}} ）}^{- β_{p}}} ．$

这个问题有一些物理参数:

动脉顺应性 $c_{一个} ＝ 1 ． 55 毫升 / 毫米汞柱$
静脉合规 $c_{v} ＝ 519 毫升 / 毫米汞柱$
外周阻力 $R ＝ 1 ． 05 （ 0 ． 84 ）毫米汞柱年代 / 毫升$
静脉流出阻力 $r ＝ 0 ． 068 毫米汞柱年代 / 毫升$
中风的体积 $V_{str} ＝ 67 ． 9 （ 77 ． 9 ）毫升$
典型平均动脉压 $P_{0} ＝ 93 毫米汞柱$
$α_{0} ＝ α_{年代} ＝ α_{p} ＝ 93 （ 121 ）毫米汞柱$
$α_{H} ＝ 0 ． 84 {证券交易委员会}^{- 2}$
$β_{0} ＝ β_{年代} ＝ β_{p} ＝ 7$
$β_{H} ＝ 1 ． 17$
$γ_{H} ＝ 0$

系统受外围压力的影响很大，外围压力从 $R ＝ 1 ． 05$ 来 $R ＝ 0 ． 84$ 开始在 $t ＝ 600$ ．因此，系统在低阶导数处具有不连续性 $t ＝ 600$ ．

恒解历史是根据物理参数定义的

$P_{一个} ＝ P_{0}, P_{v} （ t ）＝ \frac{1}{1 + \frac{R}{r}} P_{0}, H （ t ）＝ \frac{1}{{房车}_{str}} \frac{1}{1 + \frac{r}{R}} P_{0} ．$

要在MATLAB®中求解这个方程组，在调用时滞微分方程求解器之前，需要对方程、参数、延迟和历史进行编码dde23，这是指具有恒定时延的系统。您可以将所需的函数作为本地函数包含在文件的末尾(如这里所做的)，或者将它们作为单独的命名文件保存在MATLAB路径的目录中。

定义物理参数

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

p.ca = 1.55;p.cv = 519;’变为贬义词= 1.05;’变为贬义词= 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;

代码延迟

接下来，创建一个变量τ表示恒定时间延迟 $τ$ 在这些项的方程中 $P_{一个}^{τ} （ t ）＝ P_{一个} （ t - τ ）$ ．

τ= 4;

代码方程

现在，创建一个函数来编码这些方程。这个函数应该具有签名dydt = ddefun (t, y, Z, p),地点:

t为时间(自变量)。
y是解(因变量)。
Z (n, j)接近延迟 $y_{n} （ d （ j ））$ ，其中延迟 $d （ j ）$ 由分量给出j的德(t, y)．
p是可选的第四个输入，用于传入参数的值。

前三个输入由求解器自动传递给函数，变量名决定如何编写方程。参数结构p在调用求解器时传递给函数。在这种情况下，延迟表示为:

Z (: 1) $\to P_{一个} （ t - τ ）$

函数dydt = ddefun (t, y, Z, p)如果t <= 600 p.R = 1.05;其他的pr = 0.21 * exp(600-t) + 0.84;结束ylag = Z (: 1);Patau = ylag (1);Paoft = y (1);Pvoft = y (2);Hoft = y (3);dPadt = - (1 / (p.ca * p.R)) * Paoft.．.+ (1 / (p。c一个* p.R)) * Pvoft.．.+ (1/p.ca) * p.Vstr * Hoft;dPvdt = (1 / (p.v v * p.R)) * Paoft.．.- 1 / (p.cv * p.R).．.+ 1 / (p.cv * p.r)) * Pvoft;Ts = 1 / (1 + (Patau / p))^p。贝塔);Tp = 1 / (1 + (p.alphap / Paoft)^p。betap);dHdt = (p.alphaH * Ts) / (1 + p.gammaH * Tp).．.- p. h * p;dydt = [dPadt;dPvdt;dHdt];结束

注意:在示例的最后，所有函数都作为局部函数包含在内。

代码解决方案历史

接下来，创建一个向量来定义这三个组件的恒定求解历史 $P_{一个}$ , $P_{v}$ , $H$ ．解决历史就是解决时间 $t \leq t_{0}$ ．

P0 = 93;Paval = P0;Pvval = (1 / (1 + p.R/ p.R)) * P0;Hval =(1 /(’变为贬义词* p.Vstr)) *(1 /(1 +’变为贬义词/’变为贬义词))* P0;历史= [Paval;Pvval;Hval];

解决方程

使用ddeset说明不连续的存在 $t ＝ 600$ ．最后，定义积分区间 $［ t_{0} t_{f} ］$ 求解DDEdde23解算器。指定ddefun使用匿名函数传入参数结构，p．

选择= ddeset (“跳跃”, 600);Tspan = [0 1000];sol = dde23(@(t,y,Z) ddefun(t,y,Z,p)， tau, history, tspan, options);

策划解决方案

解决方案结构索尔的字段sol.x而且sol.y其中包含求解器所采取的内部时间步骤以及在这些时间点对应的解。(如果你在特定的地方需要解决方案，你可以使用德瓦尔在特定点上求解。)

绘制第三个溶液组分(心率)与时间的关系图。

情节(sol.x sol.y(3:))标题(“压力反射反馈机制的心率。”)包含(“t”) ylabel (“H (t)”）

图中包含一个axes对象。坐标轴对象的标题为心跳为反射反馈机制。包含line类型的对象。

本地函数

这里列出的是DDE解算器的本地帮助函数dde23调用来计算解。或者，您可以将这些函数作为它们自己的文件保存在MATLAB路径下的一个目录中。

函数dydt = ddefun (t, y, Z, p)%方程式被解如果t <= 600 p.R = 1.05;其他的pr = 0.21 * exp(600-t) + 0.84;结束ylag = Z (: 1);Patau = ylag (1);Paoft = y (1);Pvoft = y (2);Hoft = y (3);dPadt = - (1 / (p.ca * p.R)) * Paoft.．.+ (1 / (p。c一个* p.R)) * Pvoft.．.+ (1/p.ca) * p.Vstr * Hoft;dPvdt = (1 / (p.v v * p.R)) * Paoft.．.- 1 / (p.cv * p.R).．.+ 1 / (p.cv * p.r)) * Pvoft;Ts = 1 / (1 + (Patau / p))^p。贝塔);Tp = 1 / (1 + (p.alphap / Paoft)^p。betap);dHdt = (p.alphaH * Ts) / (1 + p.gammaH * Tp).．.- p. h * p;dydt = [dPadt;dPvdt;dHdt];结束

参考文献

[1] Ottesen J. T. <带时滞的barreflexx -反馈机构的建模>。j .数学。医学杂志。1997年第1期第36卷第41-63页。

另请参阅

ddensd|ddesd|dde23|德瓦尔