Python：如何从自定义ODE函数提取非微分方程变量的全量值？

这个问题我碰到过，其实有几种靠谱的解决办法，我给你详细说说：

方法1：求解完成后重新计算（最推荐） #

odeint返回的sol数组已经包含了每个时间点上A、B、C的数值，我们完全可以利用这些结果，事后批量计算所有时刻的Var值——这是最稳妥、最清晰的方式，还不会污染原有的微分方程函数。

先修正你代码里的小问题（比如init长度和z的索引不匹配，还有k未定义），然后实现这个方法：

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

# 先定义k的取值，你代码里没给出，这里假设为1
k = 1

def func(z, W):
    A = z[1]
    B = z[2]
    C = z[3]
    Var = k * B / C
    dAdW = A * 3
    dBdW = B * 2 / Var
    dCdW = C ** 2 * Var ** 2
    return dAdW, dBdW, dCdW

# 修正init长度，因为z用到了z[3]，所以需要4个初始值
init = [0, 0, 1, 1]
W = np.linspace(0, 100, 1000)
sol = odeint(func, init, W)

# 从sol中提取每个时间点的B和C，批量计算Var
B = sol[:, 2]
C = sol[:, 3]
Var = k * B / C

# 现在可以正常绘图了
plt.plot(W, sol[:, 0])
plt.plot(W, Var)
plt.show()

这个方法的好处：

• 原函数func保持纯净，没有额外副作用，方便后续测试和复用
• 计算逻辑直观，不会因为ODE求解器的内部迭代逻辑出现数据不匹配的问题
• 调试和维护都很简单

方法2：用可变对象实时存储中间值（适合特殊场景） #

如果必须在求解过程中就记录Var的值，可以在外部定义一个列表，每次调用func时把当前的Var追加进去。注意要用可变对象（比如列表），因为不可变对象在函数内部修改不会影响外部：

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

k = 1
# 定义一个列表来存储所有Var值
var_history = []

def func(z, W):
    global var_history
    A = z[1]
    B = z[2]
    C = z[3]
    Var = k * B / C
    var_history.append(Var)
    dAdW = A * 3
    dBdW = B * 2 / Var
    dCdW = C ** 2 * Var ** 2
    return dAdW, dBdW, dCdW

init = [0, 0, 1, 1]
W = np.linspace(0, 100, 1000)
sol = odeint(func, init, W)

# 注意：odeint可能会为了计算精度额外调用func，导致var_history长度大于W
# 所以取前len(W)个值来匹配时间点
plt.plot(W, sol[:, 0])
plt.plot(W, var_history[:len(W)])
plt.show()

⚠️ 这个方法有个坑：odeint在自适应步长求解时，会多次调用func计算中间步，导致var_history的长度比W长，需要额外处理。所以除非你有特殊需求，不然还是优先用方法1。

方法3：改用solve_ivp（更现代的ODE求解API） #

odeint是比较老旧的API，scipy现在更推荐用solve_ivp，它支持回调函数，可以精确记录每个指定时间点的Var值，不会有多余的中间数据：

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

k = 1
var_history = []
W = np.linspace(0, 100, 1000)

# solve_ivp的函数参数顺序是(时间, 状态变量)，和odeint相反
def func(W, z):
    A = z[1]
    B = z[2]
    C = z[3]
    Var = k * B / C
    dAdW = A * 3
    dBdW = B * 2 / Var
    dCdW = C ** 2 * Var ** 2
    return [dAdW, dBdW, dCdW]

# 回调函数：每个指定时间点都会调用，记录Var
def record_var(W, z):
    B = z[2]
    C = z[3]
    Var = k * B / C
    var_history.append(Var)
    return False  # 返回False表示继续求解

init = [0, 0, 1, 1]
# t_eval指定要输出的时间点，callback指定回调函数
sol = solve_ivp(func, [W[0], W[-1]], init, t_eval=W, callback=record_var)

# 现在var_history的长度和W完全一致
plt.plot(sol.t, sol.y[0])
plt.plot(sol.t, var_history)
plt.show()

这个方法适合需要在求解过程中实时处理中间变量的场景，功能更灵活。

总结一下：优先用方法1，简单可靠无副作用；如果需要实时记录，方法3的solve_ivp回调是更好的选择。

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