Handling unstable systems

[agerlach]

Unstable systems are causing extremely long quadrature run times. Seems like we can mitigate this with checking the status of the ODEProblem sol.

How should we handle the pull-back with unstable systems? Simply assign Inf to the pull back? Or assign Inf to all states in S(x) and then evaluate the observable.

[ChrisRackauckas]

How should we handle the pull-back with unstable systems? Simply assign Inf to the pull back? Or assign Inf to all states in S(x) and then evaluate the observable.

Yes, treat it as Inf. sol.retcode != :Success && sol.retcode != :Converged is a way to do this.

[aml5600]

How should we handle the pull-back with unstable systems? Simply assign Inf to the pull back? Or assign Inf to all states in S(x) and then evaluate the observable.

The latter? would allow g(x=Inf) = 0 to proceed, for example.

[agerlach]

I don't know if its so simple. What about the following

function eom(u,p,t)
  [-1.0, 1.0].*u
end

g(sol) = sol[1,end]

We have decoupled stable and unstable linear system. But g is only a function of stable state. B/c of the instability, I assume we cannot trust the value of either state and I don't know if we would even know what state is unstable. So, if I assign all the states to Inf, I would get the wrong result.

[aml5600]

I guess the question is, do we want to target the system we intended to represent via eom of eom itself? If the latter, and we consider Inf as a "member" of the state-space then the expressions above are valid.

[agerlach]

I guess the question is, do we want to target the system we intended to represent via eom of eom itself?

I'm not clear on what you are asking here.

[aml5600]

Floating point representation will send u[2] to Inf in finite t (actual system), u[2] -> Inf as t ->Inf for ODE system (probably system the user was trying to represent)

[agerlach]

Yes, but u[1]->0 as t->Inf and g() is only a function of u[1]. I interpret that you want to do g([Inf, Inf]) which isn't correct as the first state will be finite. The issue is that the integrator would have stopped early so we don't actually know what value it would be at the final time. I don't know if I have any better idea, just trying to highlight the issue.

Its seems like in the rest of the ecosystem the suggestion is to set g()=Inf: Link

The issue I see is that in the case in the link its up to the user to decide the appropriate behavior. However, here we are imbedding that decision into the expectation() call.

[aml5600]

to clarify, in the case of the two dimensions I was not saying send both to Inf, rather if one (or all) states go to Inf, send them to g as is. Wouldn't that then be up to the user?

[agerlach]

When I faced this problem, I was think it was going to actually require changes in expectation() outside of the user supplied g(). But now that I look at it more, I don't know if that is the case.

[agerlach]

Using Quadrature.jl with g(sol)=Inf leads to a bad time 😞. I was hoping the algorithms would just stop and return Inf.