我在创建一个自定义的Python类 ode_lotka，用于求解一个修正的Lotka-Volterra模型，用于描述一个群落中细菌的生长（包含一个额外的质粒转移项）。

Lotka-Volterra Model of Predator-Prey Relationship ... Lotka-Volterra（LV）模型用于细菌生长，是从生态系统中的捕食-被捕食关系模型改编的一组微分方程，用于量化群落中微生物种群之间的相互作用（竞争、促进）及其生长。它通过将指数增长与不同细菌之间的相互作用系数（竞争、共生）结合起来，预测种群数量的变化。

我已经写了函数 lotka_drug_plasmid，并使用 scipy 函数 solve_ivp 来求解它。我想在求解过程中每个时间步从这个函数中检索中间变量 plasmid_transfer。

我找到了 这个解决方案 并尝试实现，但我发现来自 solve_ivp 函数的结果与我独立输出的任何变量都不匹配。

正如 这个解决方案 所述，我尝试了：

1. 将 plasmid_transfer 与它们的求值时间一起添加到一个类字典中，然后使用 [self.plasmid_transfer_val[i] for i in model.results.t] 仅筛选被接受的时间步
2. 对每个解的时间点使用 ode_lotka.simple_plasmid(model.results.y, conjugation_rate) 事后重新计算 plasmid_transfer

这两种方法都给出稍有不同的结果，在检查输出的数值时，似乎并不正确（凭直觉来看），因为在某些时间点，plasmid_transfer 大于由 model.results.y[2] 输出的总种群数量。为了再次确认，我按照与 plasmid_transfer 相同的方式保存了函数返回值，并与输出的 model.results.y[2] 进行比较，尽管据说它们应该相同，但这些数字并不匹配。

我不确定是什么原因会导致这些结果不同，以及如何检索 plasmid_transfer 项并确保该数值确实正确。

一个简化版本的代码：

class ode_lotka():

    def __init__(self, alpha, r):
        self.alpha = alpha
        self.r = r
        self.plasmid_transfer_val = {}  # variable to retrieve intermediate term
        self.dTdt = {}  # variable to retrieve intermediate term

    @staticmethod
    def simple_plasmid(X, conjugation_rate):
        return X[1] * (conjugation_rate*(X[0] + X[2]))  # calc intermediate term

    def lotka_drug_plasmid(self, t, X, conjugation_rate, susceptibility, drug_level):
        # Lotka-Volterra equation
        dXdt = X * (self.r + np.matmul(self.alpha, X)) # calculate growth rates
        # Plasmid transfer        
        plasmid_transfer = self.simple_plasmid(X, conjugation_rate)  # calc plasmid
        self.plasmid_transfer_val.update({t: plasmid_transfer})
        # Drug susceptibility
        drug_effect = np.array(susceptibility) * drug_level
        # Define specific derivatives accounting for drugs + plasmids and return
        dDdt = dXdt[0] - drug_effect[0] * X[0]                     # donor
        dRdt = dXdt[1] - drug_effect[1] * X[1] - plasmid_transfer  # recipient 
        dTdt = dXdt[2] - drug_effect[2] * X[2] + plasmid_transfer  # transconjugant 
        dCdt = dXdt[3:] - drug_effect[3:] * X[3:]                  # community 
        self.dTdt.update({t: X[2] + dTdt})
        return np.concat([[dDdt], [dRdt], [dTdt], dCdt], axis=0)

    def solve_ode(self, t_start, t_end, x0, conjugation_rate, susceptibility, drug_level):
        self.results = solve_ivp(self.lotka_drug_plasmid,
                                 [t_start, t_end], 
                                 x0, 
                                 first_step=1, max_step=1, 
                                 args=[conjugation_rate, susceptibility, drug_level])

解决方案

你所看到的行为是预期之中的，源自 solve_ivp 在内部的工作原理，并不一定是你的模型有错。

关键点在于 solve_ivp 不会仅在被接受的时间步（results.t）对你的常微分方程进行求解。它使用自适应方法，因此你的函数在中间点（甚至在被拒绝的步长上）会被多次调用。

因此，当你像这样存储数值时：
self.plasmid_transfer_val.update({t: plasmid_transfer})
你是在收集来自所有内部求解器评估的值，而不仅仅是最终被接受的解点。这些 t 值不会与 results.t 完全匹配（由于浮点差异和被拒绝的步长），这也是为什么后续对它们进行筛选时会得到不一致或不正确的结果。

为什么你的方法会不同

1.将数据存储在ODE函数内部

• 包含中间步骤和被拒绝的步长
• t 值与 results.t 不完全对齐
• 同一时间可能会被多次求值

因此这些数据无法可靠地与最终解相匹配。

2.事后重新计算

这是正确的方法，但必须正确执行。你的函数期望一个单一的状态向量，而 results.y 是一个二维数组。

你应该这样计算：

plasmid_transfer = [
    ode_lotka.simple_plasmid(model.results.y[:, i], conjugation_rate)
    for i in range(len(model.results.t))
]

这确保了：

• 仅使用被接受的解点
• 数值与 results.y 一致
• 不会出现浮点匹配问题

关于“错误”的数值

看到 plasmid_transfer 大于 X[2] 并不一定错误。表达式

X[1] * (conjugation_rate * (X[0] + X[2]))

表示一个速率，而不是一个种群数量，因此它可以超过单个区室的容量。

另外请注意你代码中的这行：

self.dTdt.update({t: X[2] + dTdt})

与 results.y[2] 不能直接比较，因为 dTdt 是一个导数。你需要在时间上对它进行积分（即乘以一个时间步长），才能与状态值比较。

推荐的方法

最可靠的解决方案是：

求解系统后，再从 results.y 计算任何中间量。

除非你明确处理求解器内部实现，否则避免在ODE函数内部存储中间值。