在solve_ivp中缓存以前的值

Question

我想使用时间步 t=t 处的温度用于下一个时间步，因此在时间步 t=t+dt 处。我不想使用全局，因为我认为它会减慢代码速度并且不太干净。

我使用 scipy.solve_ivp 来求解质量平衡，并且函数中有一个 if 语句来检查浓度是否超过限制。当超过上限时，温度会发生变化，并且必须保持该温度，直到浓度超过下限，然后又必须变回原来的温度，依此类推。因为我在计算过程中使用函数中的温度，所以它必须使用最后一个时间步的温度。我通过将温度 T 设置为全局来解决这个问题，但我正在寻找更好的解决方案？

这是我的代码的简化版本：

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))

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

在函数中将温度 T 设置为全局变量可以使代码正常工作。我试图让它也可以通过将 T 作为参数添加到求解器中来工作，但我没有设法让它工作。

Answer 1

我最近看到了一个不同的答案，展示了如何做到这一点，但我不记得那是哪个答案（如果我再次找到它，我将使用参考更新它）。

你可以做的是利用 python 如何处理函数中的默认可变参数（这将在here讨论）。简而言之，如果您有一个可变的默认参数（例如列表），则该列表将在第一次函数调用时创建，并将在以后的所有调用中持续存在（包括更改）。例如：

def f(times_called=[0]):
    times_called[0] += 1
    print(f"This function has been called {times_called[0]} time(s)")

for i in range(10):
    f()

输出：

This function has been called 1 time(s)
This function has been called 2 time(s)
This function has been called 3 time(s)
This function has been called 4 time(s)
This function has been called 5 time(s)
This function has been called 6 time(s)
This function has been called 7 time(s)
This function has been called 8 time(s)
This function has been called 9 time(s)
This function has been called 10 time(s)

每次调用函数时，

times_called

变量都会更新，您不必处理全局变量。

所以，这里有一个与您想要的类似的示例：

from scipy.integrate import solve_ivp

def ode(t, y, y_prev=[None]):
    if t == 0:
        y_prev[0] = 0
    dydt = -y_prev[0]*y      # -y_previous*y_current
    y_prev[0] = y
    return dydt

sol = solve_ivp(ode, [0, 10], [1])

这里的主要“困难”是为

y_prev

设置一个变量初始值。您可以处理此问题的一种方法是使用

functools.partial

来填写

y_prev0

参数。

from scipy.integrate import solve_ivp
from functools import partial

def ode(t, y, y_prev0, y_prev=[None]):
    if t == 0:
        y_prev[0] = y_prev0
    dydt = -y_prev[0]*y      # -y_previous*y_current
    y_prev[0] = y
    return dydt

y_prev0 = 0
ode_partial = partial(ode, y_prev0=y_prev0)

sol = solve_ivp(ode_partial, [0, 10], [1])