I have an ODE that I can solve, but I need to also be able to solve its “inverse”. That is, my normal ode is defined as du_dx = f(x, u, p) , and I can solve that, but I also need to be able to solve dx_du = g(x, u, p) (because I’m trying to solve a boundary value problem w/ a shooting algorithm).
I’m hitting a weird boundary error when I try to do the inverse problem though, and I’m stuck on why. Any advice on what to try/check or what might help resolve this would be really helpful (and greatly appreciated).
One thing that seems to be contributing – the ODE solution should be monotone, but for some values, it seems to break and go “backward” – rather than sampling increasing alpha values, it samples a value that is below the initial value at some point and then blows up.
I’m not sure why this occurs or how to prevent it…
I can change da_ds = 0.00001 whenever ds_da = 0, but this doesn’t really help – the solver never gets past the blow up point.
You compute the flow of f as u(x). I assume everything is scalar. To compute the inverse function of u, the vector field f must not vanish. In case it does, g = 1/f will change sign through a singularity and you need to chose a solver that does not step after the singularity. Theoretically, the solution blows up at this singularity but the numerical time steps need to be adapted to not cross it.You can enforce the sign of g being constant on a trajectory using a simple implicit solver (I bet).
Not much I am afraid. CVODE_BDF() is quite good in general but I doubt it is positive (meaning it keeps the sign of the input vector). Also, around the singularity, the problem is stiff, so AutoTsit5(Rosenbrock23()) could be a first try.