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).