Update par after every time step

When I use DifferentialEquations.jl to solve a set of stiff coupled ODEs, my function is something like

function f!(du, u, par, t)
V = g(u, par)   # calculate V based on u and some other parameters
du .= ???      # some calculations to update du based on u, par and V
return nothing

du is updated based on the dependent variable u, some variables and parameters stored inside par, and also an intermediate variable V which is calculated by calling some other functions that use u and par. The problem is that calculating V=g(u,par) takes quite a long time and may require calling some external library. Since many stiff ODE solvers like the recommended FBDF requires calculating the Jacobian, it is not quite feasible to calculate the Jacobian of f! due to V. However, for my specific problem I think it might be okay that within each time step V is just treated as a variable independent of the current u, and just approximates V using the u from the previous step. Since the system should evolve to a steady state, I suspect this method is good enough for the required accuracy.

So mu question is, is it possible to just compute f! (and possibly the Jacobian) by ‘ignoring’ V=g(u,par)? Or as I suggested, is it possible that at the current time step k, calculate V=g(u_{k-1},par) where u_{k-1} is the solution at the previous time step? I’m wondering if I can have a callback to just compute V after each time step and store it in par, so that at the next time step the update of du does not require u in the current step, and so possibly the Jacobian calculation won’t take ages.

You can provide your own jacobian function and you can implement this in any way you want. For example, you could calculate a new jacobian every 10 steps or something like that.

See the options here

Use a different method (NonlinearSolve.jl) to solve to steady state?

Make V a parameter, put the g call in a DiscreteCallback, make the condition always true, and then the affect! just calls g to update the parameter V. Note that this equation won’t necessarily have a unique solution anymore since it’s time-step dependent on what V is.

The solver is already doing this :sweat: . You don’t want to mess with that because if you do it incorrectly you’ll just get divergence. This is not a good idea.