用于ODE求解的Python中的“OutputFcn”实现

Question

我正在尝试使用Python的solve_ivp来解决ODE。但是，我想基于当前解决方案与先前解决方案之间的比较动态地更改ODE的右侧。这背后的想法是我的右侧是矢量场，我想通过基于前一解决方案的方向反转右侧来确保矢量场的方向性。

其实现如下：我想检查前一个解决方案和向量字段之间右侧功能定义中的点积。如果点积为负，则右侧乘以-1。

因此，我需要访问ODE求解器的先前状态，并将其与当前迭代进行比较。在MATLAB中，有可能在解决ODE时使用“OutputFcn”。在积分器的每次迭代之后调用此函数。因此，在函数中，可以简单地将状态提取为变量，并在下一次迭代中使用它。我无法为Python找到类似的东西。

def RHS(timesnotused,x):
    out = solve_ivp(doubleGyreVar, [0,T/2, T], [x[0], x[1], 1, 0, 0, 1], rtol = 1e-10, atol=1e-10)
    output = out.y
    J = output[2:,-1].reshape(2,2)
    CG = np.matmul(J.T , J)
    lambdas, xis = np.linalg.eig(CG)
    xi_1 = xis[np.argmin(lambdas)] 
    xi_2 = xis[np.argmax(lambdas)]
    lambda_1 = np.min(lambdas)
    lambda_2 = np.max(lambdas)

    alpha = ((lambda_2-lambda_1) / (lambda_2+lambda_1))**2

    sign = 1
    if np.dot(xlast,xi_1) < 0:
        sign = -1

    return(sign*alpha*xi_1)

可以看出，我希望“xlast”成为先前的解决方案，并使用当前迭代的xi_1进行检查。不知何故，每次迭代都需要更新xlast。

Answer 1

如果要解决的微分方程是

dx/dt = sign(g(x)) * F(x)

对于一些足够平滑的函数g, F，那么你有一个不连续的右侧，所有高级数值解算器一接近这个跳跃奇点就会产生无意义。

解决这种多相系统的最明确的方法是仅向数值求解器提供连续的右侧，并通过scipy.integrate.solve_ivp中也存在的事件机制来处理相位变化。

我探讨了一种机制，用scipy的工具做类似的问题，在odeint returns wrong results for an ODE including descrete function中产生一个不连续的符号函数