How to extract intermediate calculation results from an SciPy ODE function in python?

Question

I have a bit lengthier ODE function which was simulated by using Scipy solve_ivp function. During this simulation I calculated many parameters but as the output, I am taking out put only some other parameters in return function.

But now I need to get some values of some other variables I used during the intermediate steps. So is there a way I can do that?.

Note-: When I Tried to return those parameters with the main ode outputs, it says "The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (2,) + inhomogeneous part."

Highly appreciate any valuable suggestions. Better if you can show the way by using a sample code.

Answer 1

There is no simple way of doing this. I see a few solutions though.

1. You could hack this by adding, at each function call, these secondary variables to a global list, along with their time of evaluation. Then at the end of the integration, you would need to keep only the values that correpond to the last function evaluation of each time step. Indeed, when you use methods with substeps, Jacobian evaluations or adaptive time stepping, you may have many more function calls than actual output time points.

2. Simply recompute them a posteriori for each solution time point. If you variables are simple functions of the state vector, or require simple Newton iterations, this can work quite well in my experience.

3. Modify solve_ivp do that it can accept a function which returns the time derivatives of the state vector, as well as the values of these secondary variables, and store them to put them inside the OdeSolution bunch object that is returned at the end of the integration. This is however more complex. Alternatively, solve_ivp could also call a secondary function after each accepted step, to get the corresponding secondary variables and store them.

EDIt: a cleaner and easier hack might be to use the event functionality of solve_ivp. The user-specified event functions are called only after each successful time step. This may ease the use of the first solution.