> - 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 think you should be able to use [0., Inf] in the ODE.jl solvers (or
[0, 1e308]). Of course currently, without event detection, they will
not terminate.
> 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.
Did you see DASSL.jl, it has an iterator method built in. Maybe there
you can use tspan=[0,Inf]
> 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?
This is what Alex suggested, in an example:
type ODEmeth{Name}
...
end
method_euler = ODEmeth{:euler}(...)
solver = odesolver(f, ta, ya, method_euler, deltat = 0.01)
This will dispatch on method_euler. Note, you cannot use keyword args
for dispatch.
> - 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" ?
Traits.jl might work or just use a Holy-trait...
> - 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?
I think/hope that is what ODE.jl will do in the long run as well. You can
easily go f! -> f but not f -> f! is not necessarily as performant.
> - 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
