如何在scipy.integrate.solve_ivp中缓存前步温度，避免全局变量？

问题背景 #

使用scipy.integrate.solve_ivp求解质量平衡方程时，需基于上一步的温度值计算当前时间步长。当前通过全局变量T实现温度切换逻辑：当浓度超过上限时切换为T_inside，低于下限切回初始温度T_outside，温度保持到下一次触发切换条件。但全局变量不仅可能降低运行效率，还不符合代码规范。尝试通过solve_ivp的args参数传递并更新T，未达到预期效果。

用户提供的简化代码：

import numpy as np
import time
from scipy.integrate import solve_ivp

time_check10 = True
T_outside = 20

def mass_balance(t, y):
    c = y[:len(y)//2]
    q = y[len(y)//2:]

    global T

    if t == 0:
        T = T_outside
    else:
        pass

    dc_dt = np.zeros_like(c)
    dq_dt = np.zeros_like(q)
    dc_dz = np.zeros_like(c)
    q_eq = np.zeros_like(q)

    # Boundary conditions inlet
    dc_dz[0] = (c[0] - c_feed) / dz
    # Boundary conditions outlet
    dc_dz[-1] = 0
   
    # Discretization in z-direction
    dc_dz[1:-1] = (c[1:-1] - c[:-2]) / dz
    
    # The complete mass balance equations
    dq_dt = (q_eq - q) * K * T

    dc_dt = (-u * dc_dz / epsilon) + (-((1 - epsilon) / epsilon) * rho * dq_dt)

    if q[-1] >= q_lim_outside:
        T = T_inside
    elif q[-1] <= q_lim_inside:
        T = T_outside
    else:
        pass
    
    global time_check10
    if time_check10 == True and t > 0.1*Time:
        print('10% completed after', round(time.time()-start_script,1), 'seconds')
        time_check10 = False

    return np.concatenate((dc_dt, dq_dt))

# 假设以下变量已定义
c_feed = ...
dz = ...
K = ...
u = ...
epsilon = ...
rho = ...
T_inside = ...
q_lim_outside = ...
q_lim_inside = ...
Time = ...
start_script = time.time()
t_span = (0, Time)
y_initial = ...
t_eval = ...
method = 'RK45'
dt = ...

sol = solve_ivp(mass_balance, t_span, y_initial, t_eval=t_eval, method=method, max_step=dt, dense_output=False)

优化方案：用类封装状态变量 #

全局变量的核心问题是状态管理不规范，改用类封装需要维护的状态（温度、进度标记等），既能保证状态在时间步间传递，又符合代码规范，同时避免全局变量的性能隐患。

修改后的代码 #

import numpy as np
import time
from scipy.integrate import solve_ivp

class MassBalanceSolver:
    def __init__(self, T_outside, T_inside, q_lim_outside, q_lim_inside, Time, start_script):
        self.T = T_outside  # 初始温度
        self.T_outside = T_outside
        self.T_inside = T_inside
        self.q_lim_outside = q_lim_outside
        self.q_lim_inside = q_lim_inside
        self.time_check10 = True  # 进度标记
        self.Time = Time
        self.start_script = start_script

    def mass_balance(self, t, y, c_feed, dz, K, u, epsilon, rho):
        c = y[:len(y)//2]
        q = y[len(y)//2:]

        dc_dt = np.zeros_like(c)
        dq_dt = np.zeros_like(q)
        dc_dz = np.zeros_like(c)
        q_eq = np.zeros_like(q)

        # Boundary conditions inlet
        dc_dz[0] = (c[0] - c_feed) / dz
        # Boundary conditions outlet
        dc_dz[-1] = 0
       
        # Discretization in z-direction
        dc_dz[1:-1] = (c[1:-1] - c[:-2]) / dz
        
        # 质量平衡方程（使用当前实例的温度self.T）
        dq_dt = (q_eq - q) * K * self.T
        dc_dt = (-u * dc_dz / epsilon) + (-((1 - epsilon) / epsilon) * rho * dq_dt)

        # 温度切换逻辑：直接修改实例属性，状态会保留到下一次调用
        if q[-1] >= self.q_lim_outside:
            self.T = self.T_inside
        elif q[-1] <= self.q_lim_inside:
            self.T = self.T_outside

        # 进度检查
        if self.time_check10 and t > 0.1 * self.Time:
            print('10% completed after', round(time.time() - self.start_script, 1), 'seconds')
            self.time_check10 = False

        return np.concatenate((dc_dt, dq_dt))

# 定义参数
T_outside = 20
T_inside = ...
q_lim_outside = ...
q_lim_inside = ...
Time = ...
start_script = time.time()
c_feed = ...
dz = ...
K = ...
u = ...
epsilon = ...
rho = ...
t_span = (0, Time)
y_initial = ...
t_eval = ...
method = 'RK45'
dt = ...

# 初始化求解器实例
solver = MassBalanceSolver(T_outside, T_inside, q_lim_outside, q_lim_inside, Time, start_script)

# 调用solve_ivp：用lambda包装类方法，传递固定参数
sol = solve_ivp(
    lambda t, y: solver.mass_balance(t, y, c_feed, dz, K, u, epsilon, rho),
    t_span, y_initial, t_eval=t_eval, method=method, max_step=dt, dense_output=False
)

关键说明 #

1. 状态封装：将T、time_check10等需要跨时间步维护的变量作为类的实例属性，彻底替代全局变量
2. 动态状态更新：在mass_balance方法中直接修改self.T，修改后的状态会自动保留到下一次时间步的计算中
3. 参数传递：固定参数（如c_feed、dz）通过lambda传递给类方法，避免硬编码，提升代码灵活性
4. args参数无效原因：solve_ivp的args传递的是固定初始值，每次调用微分方程函数时都会使用该初始值，无法反映中间更新后的状态，因此无法实现温度的动态切换

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