Hi Francois, Have you seen the ODE.jl package [1]? It doesn't allow you to do the things you described, but some of the github issues touch features like iterator versions of solvers or event detection. Maybe you find the discussions interesting. Comments, suggestions and/or PRs are always welcome :-)
Best, Alex. [1] https://github.com/JuliaLang/ODE.jl On Thursday, 23 April 2015 21:47:27 UTC+2, François Fayard wrote: > Hi, > > I've been looking at Julia for a long time, hoping to get some time so I can > dive in. As I am already into Fortran, C++ and Mathematica, a free dynamic > language looks like something fun to learn. Julia seems much more fun to me > than Python. > > I would like to implement an ode solver in pure Julia so we can use any type > available to Julia. What strikes me with current implementations (in > Julia/ODE, Python/Scipy, Mathematica/NDSolve) is the API that I find useless > for a program I just wrote. Let's suppose you have a Cauchy problem: y(ta)=ya > and y'(t)=f(t,y(t)). > > - In the current implementation of Julia and Python/Scipy, you need to give > all the points t1,...,tn where you want to know y at the very beginning. If > you want to solve the differential equation on [ta,tb] and you want to find a > t such that y(t)=0, you are stuck because you'll most likely do a Newton > method and your evaluation point at iteration n will depend on the value of y > at iteration n-1. Mathematica solved this problem returning an "interpolation > object" instead of some values. This is extremely useful, especially with > dense methods. > - If you want to find a t such that y(t)=0 and you don't even know an upper > bound for t (such as tb) you are stuck, even with Mathematica. > > I propose the following solution. If you want to > > solver = odesolver(f, ta, ya; method = "Euler", deltat = 0.01) > y = value(solver, t) > > solver = odesolver(t, ta, ya; method = "DenseGraggBulirschStoer", relerror = > 1.0e-6) > y = value(solver, t) > > The type of solver would depend on the method. For an Euler method, it would > just contain f, ta, ya the last t evaluated and the last value y computed. > For the DenseGraggBulirschStoer method, it would contain more information: > all the values t and the corresponding y computed, and some polynomial > coefficients for evaluation in between them. > > I have a few question to implement this idea. > - The interface makes the function odesolver not type stable as the type of > solver depends upon the method. Will it prevent static compilation when it is > released? Is there any workaround? > - For the type of y, it could by anything living in a vector space. I am > thinking of Float64, Array{Float64, 1}, Array{Float64, 2} maybe a > FixedSizedArray if such a thing exists. Is there a way to enforce that? Is > there a way to specify "any type that lives in a vector space" ? > - For the function f, I am thinking of preallocating dy_dt and call f!(dy_dt, > t, y) or call dy_dt = f(t, y). The current ODE package use the second method. > Doesn't it kill the performance with many heap allocation? Will > FixedSizeArray solve this problem? Is there a metaprogramming trick to > transform the second function into the first one? Also, the first function is > subject to aliasing between dy_dt and y. Is there a way to tell Julia there > is no aliasing? > - Some implementation might come from numerical recipes even though they also > exist in other book (not as code but as algorithm). I've seen people breaking > the copyright of numerical recipes easily. To what extend the code should be > modified so it does not break the license? Is an adaptation from one language > to another enough to say that the license does not apply ? > > Thanks
