API for ODE solvers

[acroy]

This is a summary of the discussions in #7 and #4 and intended to reflect the latest decisions about the API.

Each (adaptive) solver accepts 3 arguments (PR in #14 )

• F: the RHS of the ODE dy/dt = F(t,y), which is a function of t and y(t) and returns dy/dt::typeof(y)
• y0: initial value (no type specified -> duck typing)
• tspan: an array of points at which the solution is requested (length(tspan)>=2), the elements can be ordered ascending or descending (see ode45 should support stepping from high to low and support a vector of steps #2)

Moreover, the following keywords are supported

• norm: user-supplied norm for determining the error E (default Base.norm)
• abstol and/or reltol: an integration step is accepted if E <= abstol || E <= reltol*abs(y) (ideally we want both criteria for all solvers, done in Use keyword arguments to set the tolerance #13)
• points=:all | :specified | :final: controls the type of output according to
  • points==:all (default) output is given for each value in tspan as well as for each intermediate point the solver used.
  • points==:specified output is given only for each value in tspan.
  • points==:final output is given only for the last value in tspan, i.e. at the end of integration (speculative).
• maxstep, minstep and initstep: determine the maximal, minimal and initial integration step.

The solver returns two arrays (PR in #14)

• tout: points at which solutions were obtained (also see keyword points)
• yout: solutions at times tout; if points=:all | :specified is used, yout::Vector{typeof(y0)}. For
  points==:final the solution at the end of the integration yout::typeof(y0) is returned (speculative).

If we we decide to support non-adaptive solvers (see #9), those should have the same arguments and return values.

[tomasaschan]

About the points = :final keyword argument:

If I haven't misunderstood the discussion so far, we've assumed that points::Symbol, which would be a problem since the type of yout would then depend on the value of points, making specialization difficult. But would it not be possible to change the type of points, so that in fact one sends arguments of different type into the method? For example, if the solver is defined with a generic type parameter, we could do something like this:

abstract PointsSpecifierType

type AllPoints <: PointsSpecifierType end
type SpecifiedPoints <: PointsSpecifierType end
type FinalPoint <: PointsSpecifiertype end

function ode{T<:PointsSpecifierType}(F, tspan, y0; points::T=AllPoints, ...)
    ....
end

I'm a little unsure if this is even possible, and if so on exactly what this would mean for the syntax of the keyword argument, but I hope my idea comes across anyway. Is this approach possible? And if so, does it solve the problem of type specialization?

Another solution to get this feature without performance problems, that was suggested in #7 but not listed here, was to have an overload of ode which only gives one output parameter, so that the compiler can decide based on that:

tout, yout = ode(...) # uses current solver, with points=:specified or points=:all
                      # as allowed kwargs
y = ode(...) # no point kwarg allowed, but returns only final point (e.g. by 
             # wrapping the first one and returning yout[end])

[acroy]

@tlycken : I think the suggestion of @stevengj in #7 was based on having two functions tout, yout = odeX(F, y0, tspan; keys...) and yout = odeX(F, y0, t0, tf; keys...) with different arguments, because AFAIK dispatch on function return types is not (yet) possible (but I might be wrong).

[tomasaschan]

@acroy Ah. Well, that still leaves the first suggestion, if it's possible - to make the kwargs for points have different types, so that the types of them can be used for multiple dispatch.

[ivarne]

@tlycken It is not possible to use kwargs for dispatch. It will have to be a positional argument. One thing we can do is to make a wrapped tspan array type that we can actually use for dispatch.

y_final = odeX(F, y0, ODEFinal(t_1, t_end))

I have problems to see the value, because this is a special concept that internally will be represented like

y_final = odeX(F, y0, [t_1, t_end], points = :specified )[2][end]

And we need to present all the concepts in the last example to the user anyway.

PS: I have no opinion on the order of y0 and tspan as long as we pick something that all of Julia can use.

[tomasaschan]

@ivarne Ah, cool. It was the point about not being able to use kwargs for multiple dispatch I was forgetting. Then I think we should abandon the points = :final idea entirely; it's just not going to be worth it to get it in the API in a way that's both type-safe (i.e. performant) and simple.

(Edit: sorry for closing - a little too hot on my keyboard shortcuts...)

[acroy]

If we have to stretch too much in order to support :final it might not be worth it. Maybe we can revisit this problem at some later point ...

[ivarne]

I think that now is the right time to discuss how different features we want to have impacts the interface. I think :final might be a useful feature, but unfortunately it can not be made type stable as a keyword argument. I think type stability is a must have feature because problems that solve ODEs tend to be computationally intensive. Another suggestion is to only return the final value if tspan is a tuple, but that will also be a magical interface that is harder to discover and understand than to just appending [2][end] to the function call.

[acroy]

I still think that yout = ode!(...) could be an option, which implies changing y0. The counter arguments were that there should also be an out-of-place version and that changing y0 might not be practical.

[tomasaschan]

Do I understand correctly that the out-of-place version would then, to keep things Julian, be something like yout = ode(...) with the only difference being in not changing y0, i.e. what we're actually trying to accomplish (but can't)? If so, I think there's a risk that users try to do the latter and get confused when it doesn't work.

Overall, I think adding an example in the documentation where we use [2][end] to accomplish this will basically solve the problem, without polluting the API at all.

[ivarne]

Just want to draw attention to @JeffBezanson's comment in JuliaLang/julia@b425870 about the efficency of the t1, t2, t_cont... interface in ODE solvers and quadgk, when it is used like (odexx(f,y0,(0:0.1:100)...))

still slow as hell and still a bad idea, but at least it doesn't segfault

[stevengj]

In quadgk, you are typically only going to pass a handful of points, much smaller than the number of quadrature points, so efficiency is unlikely to be an issue in practice.

But for the ODE interface it is often desirable to use a tspan with a large number of points (e.g. for subsequent plotting or other data processing), and in this case one would be better off passing an array directly rather than splatting.

So, if we support the yfinal = ode(F, y0, t0, t1) interface at all, it might be better to support it only for two time points, returning only the final solution. Users who want more data would call ts, ys = ode(F, y0, tspan).

[tomasaschan]

Regarding performance of splattering, it's not strictly a requirement for us to use it at all - we've considered it as a way of supporting two different methods:

ode(F, y0, t0, t1, ts...; kwargs...)
ode(F, y0, tspan; kwargs...)

but there's really no reason (that I'm aware of) for us to implement the second in terms of the first, rather than the other way around. Instead of defining

ode(F, y0, tspan; kwargs...) = ode(F, y0, tspan...; kwargs...)

we could just as well (and perhaps preferrably, considering Jeff's comment) do

ode(F, y0, t0, t1, ts...; kwargs...) = ode(F, y0, [t0, t1, ts...]; kwargs...)

and avoid splattering of large arrays. The natural way to generate a big array is by using ranges (i.e. something like 0:0.01:1000) or with data which is given from outside; in both cases just passing the entire tspan without splattering is syntactically easier than adding ....

[acroy]

Basically we only need to support ode(F, y0, tspan; kwargs...). But I also wouldn't mind to have ode(F, y0, t0, t1; kwargs...)=ode(F, y0, [t0 t1]; kwargs...)[2][end] just for convenience (maybe with points=:specified to save memory).

[tomasaschan]

@acroy I'm not sure I agree we should tuck on [2][end] by default, but that's a minor detail.

It seems at least that we've reached the following conclusions here:

• The main implementation should have an API on the form ode(F, y0, tspan; kwargs...), rather than t0, t1, ts..., since calling tspan... on large ranges is not performant.
• We've decided to drop the :final option on points, in favor of type stability.

I think in a first iteration (say, for version 0.2) we could be happy without various convenience features we've discussed - neither an overload on the form ode(F, y0, t0, t1, ts...) nor a method that just returns the final point are crucial for the library to be practically usable, and since they will just wrap the main implementation internally anyway it's easy to add them later without breaking anything.

[tomasaschan]

Another feature which I'd find useful is to be able to use event functions (see #11) to terminate the solver, rather than a specified time range. In that case, I'd want to use an infinite time span and provide an event function which terminates the solver when some condition is met. (Of course, if I provide an infinite time span and a terminate condition which will never be satisfied, I have only myself to blame and Ctrl+C is my friend...)

Is this something we'd like to support? In that case, what would be a good way to specify the API for it?

[acroy]

@tlycken : I fully agree that any convenience features can be postponed and supporting something like:tfinal is mainly for convenience.

Regarding the "abort on event": In principle you can use tspan=[t0 Inf], which at the moment will run forever. The abort shouldn't be hard to support once we have an event mechanism in place.

[tomasaschan]

@acroy That's a good point - the fact that tspan = [t0 Inf] "just works" is just another case of showing that Julia is a well-implemented language. I also really like @ivarne's ideas in #27, since that would make it even simpler to grasp than the somewhat unintuitive semantics of event functions (at least the way they work in MATLAB).

[tomasaschan]

@ivarne mentioned in #27 that when solving problems for t→∞ he had some problems with y→0+ overshooting, taking him into an imaginary solution space. Would it be interesting to be able to specify some constraints on the solution?

For example, one could imagine specifying constraints = :real or :positive, and, a first approach, if (all components of) the solution y do(es) not fulfill a constraint, the step size is reduced - if hmin is reached an error is thrown. It could also be generalized to constraints = c(y) (or c(t,y)) where c is a function which returns true if the solution point is valid, and false otherwise. Using the one-argument version has the advantage that it's easy to specify e.g. constraints = Base.real and re-use many already defined functions.

[ivarne]

An important difference in Julia versus Matlab, is that you will not get imaginary results for sqrt(-1). You get a DomainError telling you to convert to complex first (eg. sqrt(complex(-1.))). I think the more obvious solution is to let the fact that F(t,y) raises an exception as a sign that the step is overshooting, and we need to take a shorter step. Maybe we want a retries or backstep kwarg to tell the solver the total number of exceptions to eat, before giving up?

[mauro3]

EDIT: as suggested by @ivarne, now a separate issue #30

I have not seen the PETSc ODE/DAE time steppers and their API
mentioned in this and related issues. I think the PETSc
framework is the nicest around and the most modern (afaik).
It has been implemented over the last few years with modern
developments in mind (they are citing papers from 2003, 2009), which
is more than can be said for most other suites of ODE solvers
which carry decades old ballast around. It supports, explicit,
implicit and IMEX methods for ODEs and DAEs all within the same
API. The API is based around formulating the problem like

F(t, u, u') = G(t, u)

For problems using implicit methods this is used as F(t, u, u') =0 for
explicit ones to u'=G(t, u) (i.e. F=u'), IMEX mixes the two at once.

Best have a look at the PETSc manual which gives a nice concise intro (Chapter 6):
http://www.mcs.anl.gov/petsc/petsc-current/docs/manual.pdf

PETSc has a liberal licence so code could probably be translated
from PETSc into Julia. The problem with that approach is that
one needs to understand the whole PETSc framework, otherwise
the code looks quite obfuscated.

Further, some cool solvers PETSc has:

• General linear methods:
  http://www.mcs.anl.gov/petsc/petsc-3.4/docs/manualpages/TS/TSGL.html#TSGL
  These methods contain Runge-Kutta and multistep
  schemes (e.g. BDF) as special cases. As far as I can tell with
  my limited ODE/DAE solver knowledge, these are the best
  currently available.
• IMEX: implicit explicit methods
  http://www.mcs.anl.gov/petsc/petsc-3.4/docs/manualpages/TS/TSARKIMEX.html
  "These methods are intended for problems with well-separated
  time scales, especially when a slow scale is strongly nonlinear
  such that it is expensive to solve with a fully implicit
  method."

Disclaimer: I have not actually used the PETSc solvers, just read
the docs and wished I could use them in Matlab...

[ivarne]

@mauro3 Thanks for your input. I would prefer if you could copy your comment into a new issue. I have tried to collect the result into https://github.com/JuliaLang/ODE.jl/blob/master/doc/api.md, and if we make changes to accomodate the methodology in the PETSc solvers, we must do it there. I think the discussion will be more focused if we have separate issues for separate concerns, and not just add random stuff at the end of general issues.

[mauro3]

@ivarne my comment above was really about the PETSc API, so I thought that it would fit into this issue. I opened a new issue for it now. Let me know if you want me to add a 'under consideration' section for PETSc API to https://github.com/JuliaLang/ODE.jl/blob/master/doc/api.md.

[tomasaschan]

@ivarne Perhaps we should close this issue now, and refer further discussions to new issues and api.md?

[ivarne]

Yes. I almost closed it 2 hours ago, but felt it would be wrong to do that before @mauro3's comments had been moved to a clean issue.