Iterator versions of ODE solvers

[ivarne]

I was thinking about #11 and just relaised that an iterator version of the ODE solvers might be a viable solution. As this is a different idea than the event function, I opened a new issue.

ode = ode_45_iter(f, y0, [0,Inf])
for t, y in ode
    if y[1] > 10
        break
    end
end
# and if we want the ODE_iter to store the intermediary solutions
# (maybe that should be left as a task for the user).
t_out, y_out = get_ode_solution_vectors(ode)

This can also be very convenient for plotting partial results and to be able to track the solvers progress.

[jiahao]

I'm very amused to see this (as well as other convergent discussions), as I've been secretly working on writing IterativeSolvers as iterators and I really like how it's turning out. See this gist for a working implementation of conjugate gradients.

[jiahao]

Also observing parallel discussions about APIs makes me think that we should work toward a closer interface for both ODE.jl and IterativeSolvers.jl, especially since many of the latter can be interpreted as Galerkin-type discretizations of the former.

[tomasaschan]

This is a fantastic idea - it solves so many problems at once. Event functions are really a quite hairy solution to a problem that this approach could possibly solve much better.

It would also be entirely possible to wrap the current implementation in an iterative approach:

function ode(F, y0, tspan)
    tout = []
    yout = []
    for t, y in ode_iter(F, y0, 0) # why would we need an end value?
        push_back(tout, t)
        push_back(yout, y)
        if t >= tspan[end]
             break
        end
  end

This has the disadvantage of not giving a solution point exactly at tspan[end], but that could be solved e.g. with a keyword argument (gandalf? ;) to the iterative solver which it shall not pass.

[ivarne]

I think the iteration interface is a good thing to explore to see if we have the flexibility we need. It seems to solve some problems, and we can probably extend the start, done, next interface for more advanced control of the iteration.

@tlycken gandalf will have to work on y values also, because it would be nice to be able to keep the y values from overshooting. I have had that problem when the value of y[1]→0 when t→∞, and negative y values brought me into imaginary space.

[tomasaschan]

@ivarne Overshooting is an entirely separate problem, to which a keyword argument with a boundary value may or may not be a good solution. (And I'll get back to it shortly.) For now, I was simply referring to the fact that when I call ode(F, y0, [t0, tfinal]) I expect to have a solution point exactly at tfinal - and without any way of controlling the step length of the iterative solver in each step, we wouldn't be able to use it as our main implementation. Thus, we need a keyword argument which, if it is given, specifies a time tf (a.k.a. gandalf) which the iterative solver must not step past; if it wants to step further than tf, the step size is reduced so that the next solution point (t,y) which is returned is exactly at tf.