WIP: Api changes

[tomasaschan]

I've started polishing the API for the existing solvers, and would like some comments both on the API decisions and on the implementations of them.

The first of these commits makes a few simple changes to the api:

• it possible to run all the solvers for scalar equations, returning scalar results

  tout, yout = ode45((t, y) -> 2t, 0:2, 0) # both tout and yout are now Vector{Float}
                                           # before, you had to wrap the last argument as [0]
                                           # and the resulting yout would have been an Nx1 Float64 array
  

• ode23 and ode45 can take tolerances as keyword arguments: ode23 takes two arguments rtol and atol, and ode45 takes one argument tol. I would rather both methods took the same arguments (preferrably by introducing relative and absolute tolerance in ode45), but I didn't want to change anything in the implementation.

• ode23 and ode45 support non-standard vector norms for determining the error, by supplying a keyword argument vecnorm which is a function (x,y) -> norm(x,y) in the vector space of the solution (see https://groups.google.com/d/msg/julia-users/PouO26Q5C-c/6UBWYBq7wAwJ for some discussion) The default is the standard norm function from Base.

The second commit introduces the the notion of a "solver object" which can be re-used. It currently only handles very basic setup (no kwargs, for example), but it might be used like this:

solver = ode(:ode23, (t,x) -> 2t, 0:2, 0)
t1, y1 = solver.solve()  # solves y' = 2t, y(0) = 0 until t = 2 using ode23
solver.solver = ode45
t2, y2 = solver.solve()  # solves the same problem using ode45, without restating the problem

In the implementation of this, I am especially uncertain on whether it's better to use symbols (:ode23) or direct references to the functions (ode23) - look at the implementation to see what I mean. It could probably be prettier =)

[tomasaschan]

See also #4

[ViralBShah]

Just saw the message on julia-users. Let's transfer this package to julialang if ok and I would also integrate it better with sundials.

I like experimenting with a solver interface. Do take a look at the original paper on ODE solvers in matlab.

[tomasaschan]

Integrating with Sundials is a good idea - as an end user, I'd probably think it natural that all the ODE stuff came in the same package, so maybe we should even merge the packages (eventually - it would be nice if this package was good enough to stand on its own legs first...).

I'll definitely take a look at the MATLAB paper.

For transferring this repo to julialang, @vtjnash is the man.

[vtjnash]

GitHub tells me I can't transfer it because I don't have admin rights on JuliaLang

[ViralBShah]

Done. Can you try again?

[vtjnash]

done. (I moved a few others too -- RPMmd, Gtk, and JuliaWebRepl)

[ViralBShah]

Thanks @vtjnash

[ViralBShah]

@stevengj @tshort Thoughts on the solver interface for ODEs?

[stevengj]

Wouldn't the usual constructor for a solver object be something like ODESolver(...)? The usual Julia convention is to make the constructor have the same name as the type.

The quadgk function uses a fastnorm function by default — fastnorm is the same as norm except for matrices, where it uses an O(N^2) norm instead of the default O(N^3) norm. Since all norms on finite-dimensional spaces differ by at most a constant factor, specifying custom norms didn't seem too critical to me. But if we do this, we should have the same interface in the ODE solvers and in the integration routines; I would just use norm=... rather than vecnorm=... (passing norm=norm works just fine in Julia). Furthermore, you should expect the supplied norm function to take only a single argument, not two arguments. (So, if the user wants the infinity norm, they could use norm=x -> norm(x,Inf). We could also change fastnorm to the infinity norm on vector and matrix types.)

I would also tend to prefer an interface like: ode45(f::Function, y0, t0, ts...; options...). Putting y0 first allows us to just use varargs for the times, rather than requiring the user to construct a vector to pass, and you can force the user to pass at least one time point (throwing a MethodError, instead of a BoundsError as you do now, if they don't). This would also be more similar to quadgk.

I would also tend to use the same keywords as quadgk (which has keyword arguments reltol=sqrt(eps), abstol=0, maxevals=10^7, order=7) where it makes sense.

[stevengj]

Also, your arguments are overtyped (a common problem in the Julia libraries). You shouldn't require them to be AbstractArray subtypes. Any type for which the requisite operations are defines (+, -, multiplication by scalars, and a norm) should be acceptable. The scalar type can be inferred from the type of t, as in quadgk.

Just don't declare a type at all for the "vectors" to be integrated; this is called duck typing.

(For example, if we get a Julia chebfun implementation up and running, then the chebfun objects will be (conceptually) infinite-dimensional vectors, and won't be an AbstractArray subtype since they won't have integer-indexed components ... but they still should be perfectly possible to integrate.)

[tomasaschan]

@stevengj Thanks for a lot of useful input!

Re. solver objects: The api for this is extremely experimental at this moment. ODESolver(...) is of course a better option. What is your opinion on the choice of Symbol vs Function variables? Is it more convenient to work with ode45 or :ode45?

Re. vector norms: The idea was based on a question in the julia-users group - I don't know if I can come up with a use case for it, but since someone wants it and it's easy to implement I didn't see why not =) norm=fastnorm seems like a sensible choice of name and default value.

Re. argument order: I think there's value here in keeping the order of the arguments the same as in MATLAB (see e.g. the MATLAB help page on ode23) to avoid confusion with users coming from that direction. I can't think of a situation where I'd ask ode45 to solve an ODE for me and not want to specify at least two time values, the start and end time (or whatever else you call your t axis) - at least not until we introduce something like MATLAB's event functions, so I don't see why ode45(f, [t0 t1], y0 is a problem. (If I want intermediates I should be able to do t0:dt:t1 instead, but that's another feature...) Still, I need two values, and I need to create some kind of vector or range from them anyway, so order doesn't seem to matter much...

Re. other keyval args: I'll gladly introduce reltol and atol in all cases where tolerance levels are needed. Perhaps there's even reason to have them in all solvers, even when they don't do anything, just to unify the API? Is there anything to use arguments such as maxevals or order here? It seems to me the order is already chosen (by choosing a solver method) and maxevals is not really relevant if we already have a lower bound on the step size.

Re. duck typing: This is good input. I'll try to work something out here, but since some of the solvers use the type argument T to allocate memory, it will require some thought to do it in a way that allows for both scalars and vectors. Might take me some time =)

[stevengj]

I think it is more convenient just to call ode45 (or oderk, see below) rather than passing symbol arguments.

I don't think we should constrain ourselves by Matlab's syntax here at the expense of internal consistency or un-Julian style. Passing multiple optional arguments as an explicit array rather than as varargs is both un-idiomatic for Julia and is inconsistent with our quadgk. Note that you can still pass arrays or ranges via the ... syntax, e.g. ode(f, y0, (0:0.1:100)...). You can require at least two times by declaring the function as ode(f::Function, y0, t0, t1, ts...; keywords...), similar to quadgk.

I agree that we might as well allow an arbitrary norm to be specified via a norm=fastnorm keyword. quadgk can be updated to do the same thing.

If an algorithm doesn't support a certain keyword, e.g. abstol, it would be better to throw an error than to silently accept it and ignore it. However, I don't see why any adaptive algorithm which supports reltol cannot also easily implement support for abstol. (I don't think we should implement non-adaptive algorithms.)

I agree that you don't need maxevals if you have a lower bound on the step size (via some keyword argument), and in the case of ODE solvers there is no good way to stop and give a meaningful result if maxeval is exceeded (unlike in quadrature)`.

Regarding order, this might be a good way of specifying the order of a Runge-Kutta rule. i.e. instead of ode45, have oderk with order=4 giving ode45, but (eventually) allow the order of the Runge-Kutta rule to be changed semi-arbitrarily (like quadgk we can eventually compute the rule on the fly and cache it, in the requested precision, including arbitrary precision support).

In quadgk I use the type of the time points to set the precision of the quadrature rule and it works well; we could do the same thing here. And you can use the similar function to allocate vectors (or whatever) based on y0.

[tomasaschan]

I do see your points about argument order - not sure I am convinced just yet, but I see your points ;) If we do make this change, it might be a good thing to just have oderk with an order keyword argument, to further emphasize that it isn't the same API as MATLAB's. I've also been thinking about implementing vode, cvode et al, so having function names after method rather than order is good.

Why don't you think we should implement non-adaptive algorithms? Because they're easy enough for each user to implement, and would take up space, or for some other reason? I think there might be merit in having things like rk4, verlet and perhaps a few others really basic ones - if nothing else then as showcase for people learning Julia. (We should at the very least make sure we have at least one symplectic method.)

I'll take a stab at these changes when I find some time, and see what comes of it.

[stevengj]

I just don't see non-adaptive algorithms as being useful in a general-purpose setting. Users know what accuracy they want, not what stepsize they want, and there doesn't seem to be much if any downside to variable-stepsize adaptive algorithms. Most users would be ill-advised to use a non-adaptive method, and I don't like the idea of putting in a method that is mostly a trap for the unwary.

Moreover, every additional algorithm that is implemented increases the code complexity and maintenance effort, so it is a bad idea to include one unless it seems genuinely useful and non-redundant.

(There seem to be plenty of adaptive symplectic integrators, so that appears to be an unrelated issue.)

For pedagogical purposes, if you want to illustrate a toy ODE integrator for people learning Julia, the right thing to do is to put it in a blog post or similar, not into the production library.

[acroy]

Hi! What is the status of this PR? I just tried the API changes and they appear to be Ok. I think it would be great to have them merged, also because finding errors is much easier when the functions have a uniform interface. I am bit undecided about the "solver object", but maybe this would be worth a separate issue/PR?

There are two minor things: you forgot to change vecnorm in ode23 and maybe x or y should be used consistently for the solution (vector) everywhere? Another thing I noticed is, that t[end] is not always tspan[end]. Not sure why this would be, but I can provide an example if you want.

Finally, FWIW, I agree with @stevengj that non-adaptive algorithms are mainly creating additional maintenance overhead. Once intermediate times are properly handled by ode* one could maybe keep ode4 as an example somewhere?

[tomasaschan]

@acroy Thanks for kicking the tires a bit! I'll take a look at your comments about vecnorm and x/y naming ASAP (hopefully this afternoon). The solver object can definitely be moved to a separate discussion - I'm not sure I want it either, but found something similar in another library (I think it was NumPy) and wanted to see if I could make something useful of it.

If you'd like to provide a test case where t[end] != tspan[end] I could take a look at why that happens and see if I can fix it.

Does anyone know of a good adaptive AND symplectic algorithm, that we could add here? I do agree that non-adaptive algorithms aren't really necessary to have in a library, but without any other symplectic algorithm we might as well keep something like a Verlet integrator here just to have an option for energy-conserving problems.

[acroy]

Here is an example, which might also be used as a test for solving a system of ODEs

d = 1
r = 1
M = [[0, d, 0]'; [-d, 0, r]'; [0, -r, 0]']
omega = sqrt(d^2 + r^2)

t,y = ODE.ode45((t,y)->M*y, 0.:4*pi/omega, [0., 0., 1.])

I get

julia> t[end]
8.0

julia> 4pi/omega
8.885765876316732

The test might be

@test maximum(abs( y[end,:]' - expm(M*t[end])*[0., 0., 1.])) < tol

Regarding an adaptive symplectic integrator, SIAM J. Sci. Comput. 18, 239 (1997) might be a start? I might also have some more recent references somewhere ... Let's open a new issue for this and maybe ask in julia-users (I vaguely remember that some celestial mechanics related question showed up at some point)?

[pao]

...non-adaptive algorithms are mainly creating additional maintenance overhead.

I'd hate to not have an implementation of RK4, in no small part because I don't feel every simulation project that uses it should have to carry it around themselves. Fixed-step solvers are still the way to go for real-time, closed-loop PIL/HIL simulation tasks (and we use it for fast-time too so the PIL/HIL work can be validated.)

[acroy]

@pao: Excuse my ignorance, but why do you need RK4 in that case? Is it because you want exactly 4 evaluations of the RHS? Otherwise, using ode45 with intermediate times specified would probably be more accurate.

[pao]

with intermediate times specified

There are no intermediate times. You're only integrating one time step into the future. Your system will have changed due to the influence of the control system for the next time step. (In Simulink, this is the ode4 fixed-step solver.)

At any rate, if people are seriously considering this it might be best as a separate issue.

[jtravs]

I also very often use fixed steppers for integrating over fixed grids - it would be useful to keep this built in, but it could be in a separate infrastructure.

[stevengj]

I think that "integrating" over fixed grids is mostly a separate problem; if there is no possibility of convergence or stepsize control, then you are really solving a discrete recurrence relation, not an ODE. (And e.g. if your fixed-grid came from some other discretization, e.g. a PDE discretization, then you want to use an integration scheme based on the order/type of your original discretization.)

[jtravs]

Yes, I agree. To provide context, I often propagate a field using a PDE using the appropriate specialist methods for that problem, but need to solve some auxiliary problem with those fields as part of the PDE RHS (for example, some coupled rate equations). To do that I use a fixed step (the step determined by the field grid, which is in turn determined by other aspects of the whole problem) ODE solver. I guess this is not the main use of the ODE package, and I can implement this elsewhere, but I wanted to point out that there are some use cases for such solvers. [Of course, if you have a better way to solve my problem, I will very happily listen to your advice!]

[acroy]

I guess one could also use RK45 with a fixed-step size and return the error estimate together with the solution vector? Let's continue the discussion about that in #9 ...

[acroy]

Regarding the API I would like to make another suggestion: if tspan is a scalar, one could assume that the ode is supposed to be solved from t=0 to t=tspan and return only the final solution as a vector. This could be handy if the solution vector is large and/or one is only interested in observables obtained from the solution.

[stevengj]

@acroy, I'm not sure I like the idea of assuming that integration ranges start at 0. Couldn't the user just use tspan=[0, T] in that case?

[stevengj]

It's generally a bad idea for the type of the return value to depend on the value (as opposed to the type) of the arguments. This kind of type-instability causes bad performance to "leak" into the surrounding code. Hence I would strongly urge against having points=:final alter the return type.

One possibility would be for the syntax to be:

• tout, yout = odeX(f, tspan, y0; kws...) returns arrays of values Vector{eltype(tspan)} and Vector{typeof(y0)}, where tspan is any iterable type (hence we don't declare the type and instead do for t in tspan or similar). The keywords can include points=:all or points=:specified (although I would prefer just allpoints=true/false since it is now a binary option).
• yfinal = odeX(f, t0, tfinal, y0; kws...) returns just the final point, with the same type as y0.

I would also tend to prefer changing the argument order to odeX(f, y0, ...), so that y0 is first, in order to leave us more flexibility for adding more of tspan-like arguments at the end, and also for similarity with quadgk. However, I don't have strong feelings about this.

[acroy]

@stevengj: Sounds very reasonable. One minor issue would be that your second odeX would have to ignore the keyword (all)points, which might be surprising. Could one use odeX! for the second behavior (which would also imply that y0 is changed to yfinal)?

[stevengj]

The second variant could throw an error if the (all)points keyword is given, so I don't think this would be a problem.

You could certainly have odeX! too, but normally in Julia this would be in addition to an out-of-place version rather than instead of it. I'm also not sure that in-place updating of y0 would be desirable in practice.

[ivarne]

Is the final functionality really useful enough that we want to complicate the interface?
Indexing finalY = y[end] will be obvious for every Julia programmer, and if you set (all)points = false it will also be memory efficient (if that matters).

[acroy]

Is the final functionality really useful enough that we want to complicate the interface?

At least we shouldn't worry too much about it now. Once the PR for the (all)points keyword is on the table, we can discuss the details and decide the fate of final.

[acroy]

Should we have a keyword minstep (or so), which allows to set the minimum integration step size? There could be an exception, ODEStepException, which is thrown when the step-size becomes too small (currently there is just a println message).

[tomasaschan]

@acroy Probably a good idea.

[FedericoV]

About a month ago, I started some work on porting to Julia some code I wrote in Python to automatically calculate the sensitivity equations analytically for a system of differential equations.

This is something I've found very useful when working on fitting ODE models to a set of timecourse data, since estimating the Jacobian using finite element methods is really numerically unstable, especially when the system is very stiff.

The downside of the approach is that the individual differential equations have to be differentiable, and I've not found a good way of identifying the actual equations within the callable function, short of writing the function in a specific way that's easy to parse.

[ivarne]

@FedericoV Take a look at https://github.com/scidom/DualNumbers.jl. It uses a different approach to numerical differentiation that does not rely on a epsilon.

PS. This seems rather off-topic for the API discussion.

[FedericoV]

Umm, good point. Should I edit it out or leave it there for time being?

[acroy]

@FedericoV : That sounds very interesting! At the moment our support for stiff is ODEs is very basic (oderosenbrock) and we would really need an adaptive solver. Sensitivity analysis would also be a great (additional) feature. As @ivarne pointed out Julia already provides interesting approaches to do numerical differentiation, which could be really useful (see also the discussion here).

I think we should continue the discussion in a separate issue ...

[acroy]

Umm, good point. Should I edit it out or leave it there for time being?

Just copy it over into a new issue.

[acroy]

It seems that MATLAB doesn't allow to set minstep... But there are two other properties we might want to support: InitialStep (default is to use a step estimated from slope of the solution at the initial time, which is currently done in ode23 but not in oderkf) and MaxStep (default: 0.1*abs(t0-tf)).

@tlycken : If you don't mind, I will copy the summary of the API discussion from one of my comments above into the description of this PR to make it easier to spot. Or should we put it somewhere else?

[ivarne]

Summarizing will be great, wherever you decide to put it.

Eventually I think we will want the reasoning behind the API in a /doc/api.md file, in a Python PEP style document. A text document describing the considerations and where and why we do things differently than Matlab will be good for people who want to use this, and others who will want to implement more algorithms.

[tomasaschan]

@acroy I've been thinking about moving this discussion somewhere else for a while - moving your summary to some "cleaner" place (e.g. a new issue) is probably a good idea. Then we can close this PR (the changes here aren't going to make it into the repo anyway) and move all future discussion there, and keep updating the head of that issue when new decisions are made.

@ivarne Documenting reasoning behind various API decisions is a great idea. I have no hard preferences on exactly how it's done, so if you think these old issues with discussions won't be "good enough" then doc/api.md is certainly a good alternative. Maybe one could even keep versions of that document, to keep track of decisions to change stuff between versions.

[ivarne]

Old issues contain way too much meta discussion when someone tries to evaluate the ODE solver capabilities in Julia, and it is hard to find exactly the information you are looking for. I use Markdown for github comments, but I don't have any preference. Readthedocs provides a nice soulution for .rst.

[acroy]

It would be great to have a document containing a description of the API and maybe even some background material on the mathematics (does that work in Markdown or rst?). For now I would say having an issue which contains the latest decisions is most convenient.

[tomasaschan]

I think this issue has played its part, so I'm closing it for housekeeping purposes.

Discussions around API decisions can continue in #20, and discussions about documentation preferably in one or more separate issues.