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
