Boundaries for solutions?

[tpoisot]

(apologies since it's not really an issue, more of a question / possible feature)

We're using ODE.jl to work on population dynamics problems, and one of the things we know about the solutions is that they should never be negative. Is there a way to specify this information to the solver? What we do currently is retry when some solutions are negative.

[pwl]

Hi,

49 has events implemented, although they are not tested that well. With the current version of this PR you can do things like

julia> ODE.ode45((t,y)->y,[1.],[0.,10],stopevent=(t,y)->(y[1]>3))

which will integrate the equation y'=y until the time reaches 10 or until the first component becomes larger then 3. In fact, this will also trigger a rootfinding algorithm, which will locate the point where y[1] crosses 3 up to some tolerance (you can change the tolerance using roottol keyword argument).

As I already mentioned, this is an experimental and untested feature, so use it at your own risk:-). We would also be grateful if you report the bugs that you find if you decide to use this PR.

[pwl]

49 also has iterator interface (which is the main point of this PR), with it you can do something like

using ODE
t0=0.; y0 = [1.]
F!(t,y,dy)=(dy[1]=y[1])
ode=ODE.ExplicitODE(t0,y0,F!)
stepper=ODE.RKStepperAdaptive{:rk45,Float64}()
opts=ODE.Options{Float64}(tstop=10.)

for (t,y) in ODE.solve(ode,stepper,opts)
    if y[1] > 3
        print((t,y))
        break
    end
end

Although it looks a bit elaborate it stops at the first instance of y[1] being larger than 3 and it doesn't try to zero in on the exact time (so you skip the rootfinding part). It should also be faster then the standard ODE routine but you have to store the results on your own (push them to a table or something like that). Still, all of this is experimental.

[pwl]

Ok, now I think I know what you meant. You are looking for a solver that would be guaranteed to preserve the positivity of your solution, and not to simply interrupt the integration if the solution becomes negative (which would already be too late). So no, we don't have integrators that can guarantee that (at least I don't know of any integrator that would preserve some arbitrary property of a solution). I think it very much depends on the equation you have: different solvers might have different effects in combination with your equation. Try an implicit solver (ode23s) for example.

There is also a method isoutofdomain in the current ODE.jl (but not in the PR I was talking about, it was temporarily disabled there), that may help you. If isoutofdomain returnes true at some point in time the step is rejected and the step size is decreased, so if you cross the zero then the integrator goes back in time and tries again with a smaller step. Perhaps redefining isoutofdomain as something like isoutofdomain(y)=y <0 | isnan(y) could solve your problems?

Sorry for the confusion.

[dpsanders]

Shouldn't this be guaranteed by the equations themselves? Are there integrators that are particularly known to satisfy this?

On 17 Jul 2016 17:34, "Timothée Poisot" notifications@github.com wrote:

(apologies since it's not really an issue, more of a question / possible feature)

We're using ODE.jl to work on population dynamics problems, and one of the things we know about the solutions is that they should never be negative. Is there a way to specify this information to the solver? What we do currently is retry when some solutions are negative.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/JuliaLang/ODE.jl/issues/109, or mute the thread https://github.com/notifications/unsubscribe-auth/AALtTptYMSuZkwJ7glL0_EqjsRWDakZ6ks5qWktsgaJpZM4JOP0Z .

[mauro3]

Note, the isoutofdomain function is only used by the Runge-Kutta solvers and that it is applied component-wise. Thus if you define isoutofdomain(y::Float64)=y <0 || isnan(y) then all components of the solution need to be non-zero. If that is not the case, maybe transforming the other components may help. If that does not work, then you need to define a custom type, see PR #107.

This begs the question: Should we make isoutofdomain(x,i) dependent on the index i as well?

@dpsanders: I don't think there are special integrators for that.

[dpsanders]

Indeed, normally in population etc models, one species being 0 (extinct) is a fixed point, which any decent integrator should not be jumping over.

[mauro3]

I don't think many (any) integrators will find the fixed points for you and do something special with them. But knowing the fixed-point it should indeed be easy to tell ODE.jl to not step over it.

[tpoisot]

@dpsanders theoretically yes, at least when there are a reduced number of populations. We work on system with 50-200 populations (so 50-200 equations), and these can give some strange results. It's rare, but it happens.

@mauro3 I like the isoutofdomain approach. Our only condition is that all y must be positive or null, so it would work for us.

[mauro3]

Cool, let us know how it goes. Probably would be good to have one such model as a test-case.

[ChrisRackauckas]

Closing since isoutofdomain seems to have fixed this.