如何在求解耦合微分方程时引入突变？

Question

我一直在尝试使用odeint和solve_ivp来求解以下三个耦合微分方程系统。(这不是确切的系统。我对它进行了修改，只是为了提高问题的可读性。)

def model(x,t):    # x is a list containing the values of u,v,w as they change with time t 
    u = x[0]
    v = x[1]
    w = x[2]
    
    dudt = u - u*v                           
    dvdt = -v + u*v - v*w**2          
    dwdt = -w + 2*v*w
    
    return [dudt, dvdt, dwdt]

这可以很好地工作。但现在我想用下面的方式修改它:当u，v或w中的任何一个低于阈值时，它被重置为零，然后让系统进化。每当这三个参数中的任何一个自动超出阈值时，都会发生这种情况。进化系统的“规则”保持不变。

我试着通过这样做来修改代码：

def model(x,t):    # x is a list containing the values of u,v,w as they change with time t 
    u = x[0]
    v = x[1]
    w = x[2]
     
    if u < u_threshold:
        x[0] = 0
    
    dudt = u - u*v                           
    dvdt = -v + u*v - vw**2          
    dwdt = -w + 2*v*w
    
    return [dudt, dvdt, dwdt]

我只为你展示了它，但是你已经明白了。这似乎不起作用。

请注意，我不能在任何变量达到阈值时都停止模拟，因为这只是一个玩具模型。稍后，我将把它推广到由数百个耦合方程组成的系统。

Answer 1

它不起作用，因为状态向量x是只读的，没有副作用将更改传输回积分器的状态向量。

即使它成功了，这也不是一个好主意

• ODE函数model大多不会在解曲线的点上被调用，然后在高阶RK方法中也不会随着时间的增加而增加。

• 这样的跳跃会使常微分方程变得高度不平滑，破坏了时间步长控制器中使用的所有假设。这样做的结果可能会有所不同，但对于最低成本的集成来说都是不利的。

所以，是的，最好的方法确实是使用solve_ivp的事件机制停止模拟，或者基于solve_ivp背后的stepper类设计您自己的时间循环。