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

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方法中也不会随着时间的增加而调用

• 这样的跳跃将使ODE高度不平滑，破坏时间步长控制器中使用的所有假设。这样做的结果可能会有所不同，但都不利于最小努力的集成

因此，是的，最好的方法确实是使用solve_ivp的事件机制停止模拟，或者根据solve_ivp后面的步进器类设计自己的时间循环