Implement ODE solvers

[ViralBShah]

A quick and easy way to get ode23 and ode45 like in MATLAB is:

http://www.netlib.org/ode/rksuite/

The state of the art is SUNDIALS, but it is no longer maintained. The lead developer went to develop video games! It is widely used in the national labs, and can give us ode15s.

https://computation.llnl.gov/casc/sundials/main.html

The long term strategy may simply be to implement these solvers in julia itself. The design of the Matlab ODE suite is described in:

http://www.mathworks.com/help/pdf_doc/otherdocs/ode_suite.pdf

[ViralBShah]

Move ode23 and ode45 from examples/ to j/ in commit 84757050b26ed549b9aee77ac7c204d9963285a2. Can close this one once we have a basic ode15s implementation.

[rgbrown]

Just a comment -- since your initial post there has been a new release of SUNDIALS (in March this year -- the first in 3 years). Worth giving it another look?

[vtjnash]

@ViralBShah can we close this, as something for an eventual package contribution (or move it to a special place for issues pertaining to desirable packages)?

[ViralBShah]

I think we can close this. This will certainly be a package.

[stevengj]

There doesn't seem to be a package for this yet, and the ODE code mentioned above seems to have disappeared from Julia.

ODE solvers are bread and butter for numerics; this really needs to be done. It is especially important because this is an area where Julia can really shine: small systems of ODEs are difficult to vectorize, so in other dynamic languages you are basically forced to drop down to C if you need it to be fast.

[aviks]

The ODE code in extras is now in a package: https://github.com/vtjnash/ODE.jl

[pao]

As there is definitely a package, going to reclose this.

[stevengj]

Sorry, got fooled by the 0.1 package listing, which is missing ODE.jl.

Are there any benchmarks (accuracy vs. number of function evaluations) of how this compares with existing ODE packages (e.g. Matlab, GSL, SUNDIALS...)? i.e. Do we have a competitive implementation or something half-baked?

[pao]

I don't think so. Viral ported some of them from Octave; I added a few others from lecture notes from a dynamics simulation course I took. I think all we know is that they don't work particularly well in general (well, except for RK4). A good testsuite with appropriate problems would be a good place to start. Rewriting the whole thing is not out of the question; Julia's changed a lot in the past year or so.

[ViralBShah]

What is the way to go forward here? In the past, I have used rksuite from netlib to get ode23 and ode45 capabilities similar to matlab, and sundials for stiff and multi-step solvers. I am not sure how actively sundials is maintained now, though. If GSL's ODE solvers are maintained, and they are considered to be good, that could be a way forward.

In any case, a good testsuite would certainly help. There are lots of ODE testsuites out there, which we could pick from, rather than develop one from scratch.

[ViralBShah]

ode23 comes from Cleve's textbook, and ode45 is ported from octave. We need a serious effort for ODE stuff - this was just some fun stuff we put along the way, and it is all half-baked.

[pao]

https://github.com/tshort/Sundials.jl exists now, too.

[ViralBShah]

Sweet!

[stevengj]

We should really have a solid (production-quality, not textbook quality) implementation of (at least) ode45 in Base, supporting arbitrary normed vector spaces (any type supporting +, -, * by Real, and norm), similar to quadgk, in addition to any external package.

An external package (especially one whose README says, "Some of these solvers are known to be produce incorrect results" — apparently not production quality!) is not really sufficient. And wrappers around C libraries (though valuable!) will never support arbitrary types (not to mention BigFloat and other general numeric types).

[vellamike]

Is there any progress on this?

[ViralBShah]

I really would like to see Sundials.jl and a good ode45 from ODE.jl in Base. No progress on this front yet, but the packages are quite usable.

[glesica]

I'm interested in putting some effort into this. I'm an MS student in CS. I'm definitely not a domain expert (although I'm not completely ignorant either), but since I can almost certainly integrate (hehe) anything I do here into my thesis work (DEM simulation in Julia) I've got some time to pour into it.

My first instinct would be to port the ODE code from SciPy (http://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html) since I've worked with it before. This group might want to go in a different direction and that would be fine too. Any suggestions or nuggets of wisdom?

[ivarne]

I would read up on the issues in https://github.com/JuliaLang/ode.jl, especially https://github.com/JuliaLang/ODE.jl/pull/7, where API discussions take place right now. The current sentiment is to remove things from Base and have a set of default packages that will be installed with Julia, so I do not think a ODE solver belongs in a shrinked core.

[ivarne]

And of course also the GSOC issue https://github.com/JuliaLang/ODE.jl/issues/18. Whether you are interested in that or not, all the advise will be perfect for someone wanting to get into the project.

[ViralBShah]

What we need is to clean up the APIs in ODE.jl, and Sundials.jl, document them, have good tests, etc. If we can do some of this, it would be nice to then have some of the pure julia ODEs in base. It would not be a whole lot of code.

[stevengj]

@ViralBShah, regarding putting this in Base, one problem is that some of the ODE.jl code is derived from Octave, as I understand it, and thus is GPLed.

Regarding unifying the API, it might also be nice to use the same API in @jiahao's GSL package.

[StefanKarpinski]

I don't really understand why we should have any ODE code in Base. This strikes me as a legacy of Matlab that doesn't really make much sense. Honestly, a lot of the linear algebra stuff that's in Base right now should probably not be.

[tknopp]

Sounds like a great candidate for a "default package" (#1906, #5155)

[jtravs]

One way around the licensing issue is to derive codes from the solvers at: http://www.unige.ch/~hairer/software.html

In particular the MATLAB ode45 is very similar to the dopri5 code at: http://www.unige.ch/~hairer/prog/nonstiff/dopri5.f and this is the same set of coefficients used in ODE.jl (the Dormand and Prince pair). These solvers are BSD licensed, are of very high quality, are theoretically explained in detail in the books by the authors: Hairer, Norsett and Wanner (1993): Solving Ordinary Differential Equations. Nonstiff Problems. 2nd edition. Springer Series in Comput. Math., vol. 8.

What is particularly nice about dopri5 is the dense output. i.e. if the ODE solver takes a natural step from t0 to t1 based on the error control, dense output is efficiently available anywhere between t0 and t1.

[stevengj]

@jtravs, I like the idea of porting dopri5.

@StefanKarpinski, I feel like an adaptive non-stiff ODE solver is one of the basic building blocks of numerical methods, along with dense linear algebra, smooth quadrature, special functions, and FFTs. Furthermore, unlike (e.g.) iterative solvers or multidimensional integration, there is pretty much one "default" way to integrate non-stiff ODEs that 99% of people use (adaptive 5th-order Runge-Kutta) and which continues to be a well-regarded method by experts. This, to me, makes ode45 a good candidate to go into Base (once an implementation stabilizes).

[StefanKarpinski]

Ok, I guess. I'm biased by what I use, which is never ODEs. For what it's worth, I also think that FFTs could really be in a package, as could smooth quadrature, and special functions. Dense linear algebra, on the other hand, I use all the time, even when I'm not doing something obviously numerical.

[acroy]

@jtravs, @stevengj : I think the implementations based on embedded RK pairs (Dormand-Prince et al) are all fairly similar. You can also find one in Numerical Recipes Chap 17.2 (including dense output), which might not be suitable license-wise. But Hairer is of course an excellent resource for all kind of solvers.

Just as an info: The current code in ODE.jl was based on an implementation for Octave (the latest version can be found here), which supported only Dormand-Prince and Fehlberg 4(5) pairs. Meanwhile our code is able to handle other pairs (Cash-Karp 4(5) and Bogacki–Shampine 2(3)) and has support for both reltol and abstol.

[stevengj]

Unfortunately, if you started with the Octave code, it is still a derived work and hence still GPL no matter how many improvements you make.

[acroy]

Yeah, I had a chat about that with @pao (who ported the code to Julia). I guess what I was trying to say is that DOPRI5 would duplicate to some extent functionality we already have (to me dense output is a very nice, but independent feature), which happens to be GPLed.

If you want to have a swiss-knife non-stiff solver in Base I would rather suggest to look into DOP853 by Hairer, which uses a different RK pair and is also recommended by NR for production work (Chapter 17.2.4).

[jtravs]

Well I've used both DOPRI5 and DOP853 extensively, and while the literature always seems to promote the latter, I usually find DOPRI5 is significantly more efficient for my problems and target tolerances. I think the second point is critical though, for high tolerances DOP853 wins.

Given that the fortran codes are almost identical between them, I don't think it is hard to have julia versions of both within the same framework - as you were previously suggesting for other RK coefficients.

BTW, dense output is not really independent - the beauty of DOPRI5 is that it has a very fortunate, high accuracy and very efficient dense output developed by Shampine (for Matlab I think), included. Other dense output routines often need further function evaluations.

BTW2: I was only noting these codes here, so that they can be considered in case a decision to rework ODE.jl was taken - I haven't got time to do it, so clearly my influence is minimal!!

[pao]

Yeah, I had a chat about that with @pao (who ported the code to Julia).

To clarify: I believe @ViralBShah did the initial port of ode45 from Octave. I refactored it afterwards. Not that it changes the licensing question, just don't want to take credit for his work.

[acroy]

@ViralBShah, @pao: Sorry! Obviously I mixed that up ...

[Nicksol]

There are disscusions developing ODE solver as GSOC project in JuliaLang/ODE.jl#18 .I have posted some ideas and looking forward for comments.

[mauro3]

This issue should probably be closed as there is https://github.com/JuliaLang/ODE.jl (and a few other independent ODE packages), where this discussion should take place.

[ivarne]

Yes I think you are right. There is no simple and obviously right algorithm for a ODE solver, so we will need packages for advanced usage anyway. Adding too much functionality in Base will also potentially be a problem for memory constrained systems, so there has been numerous suggestions that we should split up and load more things "on demand" (hopefully from a cached precompilled package).

[delta-pi-1701]

https://code.google.com/p/odepy/source/browse/odepy/irk/radau.py?r=0f4ccda975e6e95af9fff5dec36ce8f91874516d

[acroy]

@11718pauld : The discussion has partly been continued in JuliaLang/ODE.jl#18 and I think RADAU was mentioned there as well.