为什么从 SciPy solve_ivp 提取的中间变量与求解结果不一致？

2026-11-04数据小陈

应根据 model.results.y 重新计算 plasmid_transfer，使数值与输出时间点及对应状态一一对应。微分方程函数内部记录的是求解器调用时的临时结果，时间相同并不代表它对应最终输出状态。

plasmid_transfer 是瞬时转移速率，model.results.y[2] 是种群数量，不能通过比较二者判断求解是否正确。X[2] + dTdt 也不是求解器计算出的下一时刻种群数量，同样不能用于这一判断。

为什么函数内部记录的值可能不同？

solve_ivp 接收的函数定义是：

fun(t, X) -> dXdt

它返回状态对时间的导数。求解器为了推进一步，可能在多个内部阶段调用这个函数；误差过大时，也可能放弃一次试探步，缩小步长后重新计算。

因此，函数内部的 (t, X) 不一定是最终出现在 results.t 和 results.y 中的那一组时间与状态。

下面的记录方式还有一个问题：

self.plasmid_transfer_val[t] = plasmid_transfer

同一时间可能发生多次函数调用，状态也可能不同，而字典只保留最后一次写入的值。随后按 results.t 查询时，只核对了时间，没有核对计算时使用的状态。浮点数时间作为字典键，也不适合作为通用的结果对齐方式。

某些调用恰好对应最终输出状态时，记录值可能一致，但不能依赖这种实现细节。

原代码中的设置：

first_step=1,
max_step=1,

分别指定首次尝试的步长和允许的最大步长，不能保证每个接受步都等于 1。

根据输出解重新计算 plasmid_transfer

原来的 simple_plasmid 可以直接处理 results.y：

import numpy as np

sol = model.results
if not sol.success:
    raise RuntimeError(sol.message)

t = sol.t
X = sol.y

plasmid_transfer = ode_lotka.simple_plasmid(
    X, conjugation_rate
)

# 每一行依次为：时间、transconjugant 数量、瞬时质粒转移速率
output = np.column_stack((t, X[2], plasmid_transfer))

sol.y 的形状为 (状态变量数, 输出时间点数)。X[0]、X[1] 和 X[2] 分别包含对应种群在全部输出时间点上的数值。原公式会按元素计算：

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

所得 plasmid_transfer[j] 对应 sol.t[j] 和 sol.y[:, j]。这一做法适用于当前代码中的情况，即 conjugation_rate 为固定参数，转移项只依赖状态和该参数。

如果希望在指定时间点获得结果，可以在求解时传入 t_eval，然后对返回的 sol.y 使用相同公式。t_eval 指定输出时间，不会强制求解器以这些时间间隔进行内部步进。

为什么 X[2] + dTdt 不等于求解结果？

这里的 dTdt 表示：

[ \frac{dT}{dt} ]

它是变化速率，不是一次函数调用产生的种群增量，因此不能与数量直接相加。

即使使用步长 h 写成：

[ T(t+h) \approx T(t) + h\frac{dT}{dt}(t) ]

这也只是显式 Euler 方法的近似更新。原代码没有指定 method，默认使用 RK45。RK45 综合多个阶段的导数进行更新，不能用单次函数调用中的 X[2] + dTdt 复现其结果。

如果需要查看输出状态对应的导数，可以直接根据 sol.y 重新计算：

drug_effect = np.asarray(susceptibility) * drug_level

growth = X * (
    np.asarray(model.r)[:, None]
    + np.asarray(model.alpha) @ X
)

dTdt_at_output = (
    growth[2]
    - drug_effect[2] * X[2]
    + plasmid_transfer
)

这段代码假定 model.r 和 susceptibility 都是一维数组，长度与状态变量数一致。dTdt_at_output 表示各输出时间点的导数，仍不能与 X[2] 当作同一物理量比较。

plasmid_transfer 大于种群数量是否异常？

仅凭数值大小无法判断，因为二者单位不同。

从代码中它直接参与 dRdt、dTdt 的方式看，plasmid_transfer 应具有“种群数量／时间”的单位，X[2] 则具有“种群数量”的单位。如果状态使用的是浓度，对应单位就是“浓度／时间”和“浓度”。conjugation_rate 的单位也必须与此保持一致。

转移项表示由 recipient 转为 transconjugant 的速率。当前 transconjugant 数量较少时，转入速率的数值高于该种群当前数量，并不矛盾。

若想知道某段时间内累计发生了多少转移，需要计算：

[ C(t)=\int_{t_0}^{t}\mathrm{plasmid_transfer}(\tau),d\tau ]

可以对输出速率做数值积分，或者增加一个满足 dC/dt = plasmid_transfer、初始值为 0 的累计状态。模型还包含生长和药物作用，因此累计转移量通常不等于 transconjugant 种群的净增量。

如何进一步检查结果？

建议先移除微分方程函数中记录中间值的字典写入，让该函数只负责返回导数，再统一根据输出解计算需要展示的变量。

检查时可以关注以下几点：

确认 sol.success，并在失败时查看 sol.message。
核对种群顺序、参数单位，以及 alpha、r、susceptibility 的维度。
当前转移项在 recipient 方程中为负、在 transconjugant 方程中为正。这两个方程相加时，转移项应相互抵消；生长和药物项仍会影响总量。
若关心数值精度，可以收紧 rtol、atol，在相同输出时间点比较解和重新计算的转移速率。没有具体参数、初始值和结果数据，目前无法判断误差大小或参数是否合理。

根据 sol.y 重新计算能保证中间变量与输出状态一致，但不能单凭这一点判断模型、参数或求解精度是否正确。