如何用odeint求解含时变函数elastance(t)的耦合常微分方程组？

问题描述 #

尝试用scipy.integrate.odeint求解4个耦合常微分方程组，但包含时变函数elastance(t)的ODE求解结果不符合预期，单独绘制elastance(t)图像正常。

原代码 #

from scipy.integrate import odeint
import numpy as np; import matplotlib.pyplot as plt

def elastance(t):
    t = np.array(t)
    Senzaki_table_4=np.array([[28.38975, np.pi/2],[37.58583, .08367674],[21.02345, -1.486758], [7.665592, 2.865675], [4.809436, .1677238], [4.181973, 4.630239], [1.940692, 3.088379], [.5870049, -.3053668], [1.181256, 4.410703], [.84039, 3.181538], [.02259011, 1.242886], [.3071458, 4.156753], [.3226207, 2.946186]])

    hrt_rate = 95
    E_es = 2
    period = 60/hrt_rate
    t_sin = t/period*2*np.pi
    E1 = 1e-2*np.ones(t.shape)[:,None]@Senzaki_table_4[:,0].transpose()[None,:]
    E2 = np.sin(t_sin[:,None] @ np.arange(0,13)[None,:] + np.ones(t.shape)[:,None] @ Senzaki_table_4[:,1][None,:])
    
    #E00 = E1*E2
    E = E_es * np.sum(E1 * E2, 1)

    return E

def opencircuit(x, t):
    # constants
    pressure_conversion = 1333.2237
    L_la = 0.00005
    v_lv0 = 2
    v_d0 = 205+370+401
    R_c = 90/pressure_conversion
    R_d = 1200/pressure_conversion
    C_d = 8e-4*pressure_conversion
    p_la = 12
    initial_v_lv = 128.35607011204
    initial_v_d = 1053.96962417265
    L_la=0.00005
    L_lv=0.000416

    initial_t = 0
    num_cycles = 2.75
    hrt_rate = 95
    final_t = 60 / hrt_rate * num_cycles

    # assgin each ODE to a vector element
    Q_m = x[0]; v_lv = x[1]; Q_a = x[2]; v_d = x[3]

    # define each ODE
    dQ_m_dt = (p_la - elastance([t])*(v_lv- v_lv0))/L_la
    dv_lv_dt = Q_m - Q_a
    dQ_a_dt = 1/L_lv*(elastance([t])*(v_lv-v_lv0) - R_c*Q_a - (v_d-v_d0)/C_d)
    dv_d_dt  = Q_a - (v_d-v_d0)/R_d/C_d

    return [dQ_m_dt, dv_lv_dt, dQ_a_dt, dv_d_dt]

# initial conditions
x0 = [0, 2, 0, 205+370+401]

# test the defined odes
# print(opencircuit(x=x0, t=0))

# declare a time vector (time window)
t = np.linspace(0, 1, 1001)
x = odeint(opencircuit, x0, t)

Q_m = x[:, 0]
v_lv = x[:, 1]
Q_a = x[:, 2]
v_d = x[:, 3]

# plot the results
plt.plot(t, elastance(t))
plt.plot(t, Q_m); plt.plot(t, v_lv)
plt.plot(t, Q_a); plt.plot(t, v_d)
plt.show()

问题分析与修正方案 #

核心问题 #

1. elastance函数输入兼容性不足：原函数强制要求输入为数组，但odeint调用opencircuit时传入的t是标量，每次调用elastance([t])会返回长度为1的数组，导致ODE计算中出现不必要的数组维度运算，可能引发数值异常。
2. 冗余计算：在opencircuit中多次重复调用elastance([t])，浪费计算资源。
3. 时间范围不匹配：原代码设定模拟2.75个心动周期，但时间向量仅覆盖0到1秒，未完成完整模拟周期。

修正步骤 #

1. 优化elastance函数，兼容标量与数组输入

修改函数，使其自动处理标量/数组输入，避免手动包装数组：

def elastance(t):
    t_arr = np.asarray(t)
    Senzaki_table_4=np.array([[28.38975, np.pi/2],[37.58583, .08367674],[21.02345, -1.486758], [7.665592, 2.865675], [4.809436, .1677238], [4.181973, 4.630239], [1.940692, 3.088379], [.5870049, -.3053668], [1.181256, 4.410703], [.84039, 3.181538], [.02259011, 1.242886], [.3071458, 4.156753], [.3226207, 2.946186]])

    hrt_rate = 95
    E_es = 2
    period = 60/hrt_rate
    t_sin = t_arr/period*2*np.pi
    E1 = 1e-2*np.ones(t_arr.shape)[:,None]@Senzaki_table_4[:,0].transpose()[None,:]
    E2 = np.sin(t_sin[:,None] @ np.arange(0,13)[None,:] + np.ones(t_arr.shape)[:,None] @ Senzaki_table_4[:,1][None,:])
    
    E = E_es * np.sum(E1 * E2, 1)
    
    # 输入为标量时返回标量，避免数组维度问题
    if np.isscalar(t):
        return E.item()
    return E

2. 简化opencircuit中的时变函数调用

在ODE函数中仅计算一次elastance(t)，并复用计算结果：

def opencircuit(x, t):
    # constants
    pressure_conversion = 1333.2237
    L_la = 0.00005
    v_lv0 = 2
    v_d0 = 205+370+401
    R_c = 90/pressure_conversion
    R_d = 1200/pressure_conversion
    C_d = 8e-4*pressure_conversion
    p_la = 12
    L_lv=0.000416

    # 分配状态变量
    Q_m = x[0]; v_lv = x[1]; Q_a = x[2]; v_d = x[3]

    # 仅计算一次当前时刻的弹性值
    E_t = elastance(t)
    p_lv = E_t * (v_lv - v_lv0)

    # 定义各ODE
    dQ_m_dt = (p_la - p_lv)/L_la
    dv_lv_dt = Q_m - Q_a
    dQ_a_dt = (p_lv - R_c*Q_a - (v_d - v_d0)/C_d)/L_lv
    dv_d_dt  = Q_a - (v_d - v_d0)/(R_d * C_d)

    return [dQ_m_dt, dv_lv_dt, dQ_a_dt, dv_d_dt]

3. 调整时间向量以覆盖完整模拟周期

根据设定的2.75个心动周期，计算正确的终止时间：

hrt_rate = 95
num_cycles = 2.75
final_t = 60 / hrt_rate * num_cycles
# 生成1000Hz采样率的时间向量
t = np.linspace(0, final_t, int(final_t * 1000))

完整修正代码 #

from scipy.integrate import odeint
import numpy as np; import matplotlib.pyplot as plt

def elastance(t):
    t_arr = np.asarray(t)
    Senzaki_table_4=np.array([[28.38975, np.pi/2],[37.58583, .08367674],[21.02345, -1.486758], [7.665592, 2.865675], [4.809436, .1677238], [4.181973, 4.630239], [1.940692, 3.088379], [.5870049, -.3053668], [1.181256, 4.410703], [.84039, 3.181538], [.02259011, 1.242886], [.3071458, 4.156753], [.3226207, 2.946186]])

    hrt_rate = 95
    E_es = 2
    period = 60/hrt_rate
    t_sin = t_arr/period*2*np.pi
    E1 = 1e-2*np.ones(t_arr.shape)[:,None]@Senzaki_table_4[:,0].transpose()[None,:]
    E2 = np.sin(t_sin[:,None] @ np.arange(0,13)[None,:] + np.ones(t_arr.shape)[:,None] @ Senzaki_table_4[:,1][None,:])
    
    E = E_es * np.sum(E1 * E2, 1)
    
    if np.isscalar(t):
        return E.item()
    return E

def opencircuit(x, t):
    pressure_conversion = 1333.2237
    L_la = 0.00005
    v_lv0 = 2
    v_d0 = 205+370+401
    R_c = 90/pressure_conversion
    R_d = 1200/pressure_conversion
    C_d = 8e-4*pressure_conversion
    p_la = 12
    L_lv=0.000416

    Q_m = x[0]; v_lv = x[1]; Q_a = x[2]; v_d = x[3]

    E_t = elastance(t)
    p_lv = E_t * (v_lv - v_lv0)

    dQ_m_dt = (p_la - p_lv)/L_la
    dv_lv_dt = Q_m - Q_a
    dQ_a_dt = (p_lv - R_c*Q_a - (v_d - v_d0)/C_d)/L_lv
    dv_d_dt  = Q_a - (v_d - v_d0)/(R_d * C_d)

    return [dQ_m_dt, dv_lv_dt, dQ_a_dt, dv_d_dt]

# 初始条件
x0 = [0, 2, 0, 205+370+401]

# 设置模拟时间范围
hrt_rate = 95
num_cycles = 2.75
final_t = 60 / hrt_rate * num_cycles
t = np.linspace(0, final_t, int(final_t * 1000))

# 求解ODE
x = odeint(opencircuit, x0, t)

Q_m = x[:, 0]
v_lv = x[:, 1]
Q_a = x[:, 2]
v_d = x[:, 3]

# 绘图
plt.figure(figsize=(12,8))
plt.subplot(2,1,1)
plt.plot(t, elastance(t), label='Elastance(t)')
plt.legend()
plt.xlabel('Time (s)')
plt.ylabel('Elastance')

plt.subplot(2,1,2)
plt.plot(t, Q_m, label='Q_m')
plt.plot(t, v_lv, label='v_lv')
plt.plot(t, Q_a, label='Q_a')
plt.plot(t, v_d, label='v_d')
plt.legend()
plt.xlabel('Time (s)')
plt.ylabel('State Variables')

plt.tight_layout()
plt.show()

说明 #

• 修正后的elastance函数可直接处理标量t，彻底解决数组维度不匹配问题。
• 减少冗余计算，提升求解效率。
• 时间向量覆盖完整的2.75个心动周期，符合初始模拟设定。

内容的提问来源于stack exchange，提问作者u_binyong