Python 科学：用条件中断微分方程求解答案

【问题标题】：Python scientific: interrupt differential equation solving with a conditionPython 科学：用条件中断微分方程求解
【发布时间】：2010-11-19 17:04:37
【问题描述】：

我目前正在python下求解一个微分方程组，使用odeint模拟场中的带电粒子（来源来自这个package）：

time = np.linspace(0, 5, 1000)

def sm(x, t): 
    return np.array([x[1], eta*Ez0(x[0])])

traj = odeint(sm,[0,1.], time)

它工作正常，但我想在 x[0]

def sm1(x, t): 
    if x[0] < 0:
        return np.array([0, 0]) 
    else:
        return np.array([x[1], eta*Ez0(x[0])])

traj = odeint(sm1,[0,1.],time)

但我想有更好的解决方案。我找到了this，但在我看来，它修复了步数，这令人遗憾。任何建议表示赞赏。

【问题讨论】：

您的解决方案对我来说似乎很合理。你觉得它有什么问题？
您会收到警告：lsoda-- 在当前 t (=r1), mxstep (=i1) 在到达 tout 之前在此调用上采取的步骤在上述消息中，I1 = 500 在上述消息中，R1 = 0.4223349048304E+00 在此调用上完成的工作过多（可能是错误的 Dfun 类型）。运行 full_output = 1 以获得定量信息。

标签： python numpy scipy

【解决方案1】：

如果您编写 odeint 函数的自定义扩展，您可以让您的函数在完成时引发特定异常。在 Python 中执行它可能会使其速度大大降低，但我认为您使用 C 或 Cython 编写相同的东西。请注意，我没有测试以下内容。

class ThatsEnoughOfThat(Exception):
    pass

def custom_odeint(func, y0, t): # + whatever parameters you need
    for timestep in t:
        try:
            # Do stuff. Call odeint/other scipy functions?
        except ThatsEnoughOfThat:
            break
    return completedstuff

def sm2(x, t):
    if x[0] < 0:
       raise ThatsEnoughOfThat
    return np.array([x[1], eta*Ez0(x[0])])

【讨论】：