MATLAB 數學應用 微分方程 時滞微分方程 具有不連續性的心血管模型DDE

本文講述了如何使用 dde23 對具有不連續導數的心血管模型求解。

方程組為：

P ˙ a ( t ) = − 1 c a R P a ( t ) + 1 c a R P v ( t ) + 1 c a V s t r ( P a τ ( t ) ) H ( t ) \dot{P}_a(t) = -\dfrac{1}{c_aR}P_a(t) + \dfrac{1}{c_aR}P_v(t) + \dfrac{1}{c_a}V_{str}(P_{a}^τ(t))H(t) P˙a​(t)=−ca​R1​Pa​(t)+ca​R1​Pv​(t)+ca​1​Vstr​(Paτ​(t))H(t)

P ˙ v ( t ) = − 1 c v R P a ( t ) − ( 1 c v R P v ( t ) ) \dot{P}_v(t) = -\dfrac{1}{c_vR}P_a(t) - (\dfrac{1}{c_vR}P_v(t)) P˙v​(t)=−cv​R1​Pa​(t)−(cv​R1​Pv​(t))

H ˙ ( t ) = α H T s 1 + γ H T p − β H T p . \dot{H}(t) = \dfrac{α_HT_s}{1 + γ_HT_p} - β_HT_p. H˙(t)=1+γH​Tp​αH​Ts​​−βH​Tp​.

T s T_s Ts​ 和 T p T_p Tp​ 的項分别是同一方程在有時滞和沒有時滞狀态下的變體。 P τ a P_τ^a Pτa​ 和 P a P_a Pa​ 分别代表在有時滞和沒有時滞狀态下的平均動脈壓。

此問題有許多實體參數：

• 動脈順應性 c a = 1.55   m l / m m H g c_a = 1.55 ml/mmHg ca​=1.55 ml/mmHg
• 靜脈順應性 c v = 519   m l / m m H g c_v = 519 ml/mmHg cv​=519 ml/mmHg
• 外周阻力 R = 1.05 ( 0.84 ) m m H g   s / m l R = 1.05(0.84)mmHg s/ml R=1.05(0.84)mmHg s/ml
• 靜脈流出阻力 r = 0.068   m m H g   s / m l r = 0.068 mmHg s/ml r=0.068 mmHg s/ml
• 心博量 V s t r = 67.9 ( 77.9 )   m l V_{str} = 67.9(77.9) ml Vstr​=67.9(77.9) ml
• 典型平均動脈壓 P 0 = 93   m m H g P_0 = 93 mmHg P0​=93 mmHg
• α 0 = α s = α p = 93 ( 121 )   m m H g α_0 = α_s = α_p = 93(121) mmHg α0​=αs​=αp​=93(121) mmHg
• α H = 0.84   s e c − 2 α_H = 0.84 sec^{−2} αH​=0.84 sec−2
• β 0 = β s = β p = 7 β_0 = β_s = β_p = 7 β0​=βs​=βp​=7
• β H = 1.17 β_H = 1.17 βH​=1.17
• γ H = 0 γ_H = 0 γH​=0

該方程組受外周壓的巨大影響，外周壓會從 R=1.05 急劇減少到 R=0.84，從 t=600 處開始。是以，該方程組在 t=600 處的低階導數具有不連續性。

常曆史解由以下實體參數定義

P a = P 0 , P_a = P_0, Pa​=P0​,

P v ( t ) = 1 1 + R r P 0 , P_v(t) = \dfrac{1}{1 + \dfrac{R}{r}}P_0, Pv​(t)=1+rR​1​P0​,

H ( t ) = 1 R V s t r 1 1 + r R P 0 . H(t) = \dfrac{1}{RV{str}}\dfrac{1}{1+\dfrac{r}{R}}P_0. H(t)=RVstr1​1+Rr​1​P0​.

要在 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 來表示項 P τ a ( t ) = P a ( t − τ ) P_τ^a(t) = P_a(t−τ) Pτa​(t)=Pa​(t−τ) 的方程中的常時滞 τ。

編寫方程代碼

現在，建立一個函數來編寫方程的代碼。此函數應具有簽名 dydt = ddefun(t,y,Z,p)，其中：

• t 是時間（自變量）。
• y 是解（因變量）。
• Z(n,j) 對時滞 y n = ( d ( j ) ) y_n = (d(j)) yn​=(d(j)) 求近似值，其中時滞 d(j) 由 dely(t,y) 的分量 j 給出。
• p 是可選的第四個輸入，用于傳入參數值。

求解器自動将前三個輸入傳遞給函數，變量名稱決定如何編寫方程代碼。調用求解器時，參數結構體 p 将傳遞給函數。在本例中，時滞表示為：

• Z ( : , 1 )   → P a ( t − τ ) Z(:,1) →P_a(t−τ) Z(:,1) →Pa​(t−τ)

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 
           

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

編寫曆史解代碼

接下來，建立一個向量來定義三個分量 P a P_a Pa​、 P v P_v Pv​ 和 H 的常曆史解。曆史解是時間 t ≤ t 0 t ≤ t_0 t≤t0​ 的解。

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 來指定在 t=600 處存在不連續性。最後，定義積分區間 [ t 0   t f ] [t_0 t_f] [t0​ tf​] 并使用 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