scipy、odeint与质量守恒

Question

我在处理以下这个系统（常微分方程）时遇到了一些困难（k的值不是常量）：

def my_diff(y,t,k):
    f  = np.zeros(4)
    f[0] = - k[0]*y[0] - k[1]*y[0] - k[2]*y[0]**2
    f[1]= k[0]*y[0]
    f[2] = k[1]*y[0]
    f[3] = k[2]*y[0]**2
    return f

这里有三个化学反应正在进行：

f[0] .. concentration of raw-material
f[1] .. concentration of product 1
f[2] .. concentration of product 2
f[3] .. concentration of product 3

如果我解决这个系统，所有的东西都运转得很好，质量是守恒的。但是如果我在循环中运行这个模拟，质量就不再守恒了，而且我的错误随着循环次数的增加而线性上升。

我做的事情（代码片段）：

# solve the ODE
sol = integrate.odeint(my_diff,y,t,(k,))
# update initial conditions and solve again
y = [ sol.T[0][-1] + new_pulse,
      sol.T[1][-1] , sol.T[2][-1] , sol.T[3][-1]]

不幸的是，系统中的总质量在增加。我反复检查我的代码，但找不到任何错误。我尝试使用IDA求解器，并通过添加一个代数项来限制总浓度，但我在定义一致的初始条件时遇到了困难。

你觉得在这种模型和常微分方程求解器下，错误会很大（在20个循环后高达50%）吗？还是我应该继续寻找bug？

Answer 1

嗯，发帖之前先去睡觉吧……

我找到了错误。我计算体积增加的方法错了。这个常微分方程系统运行得很好。其实不需要质量平衡。现在没有错误传播的问题了。

Answer 2

你说“我的错误随着循环线性增加”。实际上，你在每个循环中都在以线性的方式增加质量，因为你在每个循环中都在给你的raw_material变量增加相同的数量new_pulse，而raw_material在你系统的总质量中是一个线性项。

我猜问题出在动态方程上，而不是代码里。更准确地说，你可能需要在其他地方添加你的new_pulse变量。