Styjun

I am trying to solve an ODE using Python's solve_ivp. However, I want to change the right hand side of my ODE dynamically based on a comparison between the current solution and previous solution. The idea behind this is that my right hand side is a vector field, and I want to ensure directionality of the vector field by reversing the right hand side based on the direction of the previous solution.

The implementation for this is as follows: I want to check the dot product in the right hand side function definition between the previous solution and the vector field. If the dot product is negative the right hand side is multiplied by -1.

I therefore need to access the previous state of the ODE solver and use it in comparison with the current iteration. In MATLAB there is the possibility of using "OutputFcn" while solving an ODE. This function is called after every iteration of the integrator. In the function it is therefore possible to simply extract the state as a variable and use it in the next iteration. I have not been able to find something similar for 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)

As can be seen I want "xlast" to be the previous solution, and check it with xi_1 of the current iteration. Somehow xlast needs to be updated every iteration.

If the differential equation you want to solve is

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

for some sufficiently smooth functions g, F, then you have a discontinuous right side where all advanced numerical solvers will produce nonsense as soon as they approach this jump singularity.

The clearest method to solve such a multi-phase system is to present only continuous right sides to the numerical solver and to handle the phase change via the event mechanism that is also present in scipy.integrate.solve_ivp.

I explored a mechanism to do that with the tools of scipy in a similar problem with a sign function producing a discontinuity in odeint returns wrong results for an ODE including descrete function