Hi all, I hope this question finds you well.
I have two state vectors X1 and X2 that comprise my total state X = [X1; X2].
X1 evolves according to dX1 = f(X1,X2) and X2 evolves according to dX2 = g(X2) - h(X1,X2,dX1) (note dependence on dX1). I have a specific timestepping scheme I’d like to use for X2 where I treat g explicitly and h implicitly, and h is linear in X2 so I have a function for X2(t+dt) ready to go in terms of X2(t).
But for X1, the function f(X1,X2) is fairly complex so I like using the Julia solvers to adaptively timestep, do local error analysis etc.
I’m currently solving the system by passing dX1 = f(X1,X2) and dX2 = 0 to a julia solver, and then using a DiscreteCallback to timestep X2 manually after every timestep. However, I assume the Julia solvers will find this approach weird (since f is returning 0 for dX2) and will result in poor solving.
I’m running a test with dtmax = 1e-5 to try and ensure that X2 is being solved well, but now I’m losing most of the benefits of the lovely prewritten adaptive solvers that are so good for X1. That might work just fine, but the current ETA is about 7 hours so I thought I’d see if there’s an obvious better way of handling this situation.
Before trying this, I was just passing the rhs function dX = [dX1; dX2] to julia solvers and using a stiff solver, but I found myself running into trouble when moving certain parameters too high. E.g. julia solvers needing increasingly small (< 10^-16) timesteps, I believe for X2. I would see X2 occasionally go negative which it shouldn’t, but the IMEX scheme I’ve implemented for X2(t+1) is positive by definition (is this iterative implicit solvers vs rearranged for X2(t+1)?).