ODE.jl icon indicating copy to clipboard operation
ODE.jl copied to clipboard

ode45 hangs on backwards integration

Open crbinz opened this issue 8 years ago • 6 comments

Hello,

I came across this issue (recreated here using example code):

julia> using ODE

julia> function f(t, y)
           # Extract the components of the y vector
           (x, v) = y

           # Our system of differential equations
           x_prime = v
           v_prime = -x

           # Return the derivatives as a vector
           [x_prime; v_prime]
       end;

julia> # Initial condtions -- x(0) and x'(0)
       const start = [0.0; 0.1]
WARNING: redefining constant start
2-element Array{Float64,1}:
 0.0
 0.1

julia> time = 0.:-.1:-4pi;

julia> t, y = ode45(f, start, time);

This hangs for a long time (I've waited up to a minute before killing it).

crbinz avatar Jun 24 '16 19:06 crbinz

I belive I have the same problem. ode45 seems to go into an infitine loop under some conditions.

For example:

function f(t,v,p)
    (x,y) = v
    dx = p[1]*x - x*y
    dy = p[2]*x*y - y
    return [dx;dy]
end

p = [4.82,-2.43]
init_vals = [0.5,0.5]
times = 1.0:1.0:10.0
g = (t,v) -> func(t,v,params)
(t,d) = ode45(g,init_vals,time_points,points=:specified)

I also waited longer than a minute before killing it.

galbarel avatar Jul 01 '16 14:07 galbarel

Maybe this has to do with the 1.01 here: https://github.com/JuliaLang/ODE.jl/blob/a0ef97d6249e62a089370787850ae697e3d94962/src/runge_kutta.jl#L341 . Which I think is the wrong thing to do.

Could you check if setting it to 1 works?

mauro3 avatar Jul 01 '16 14:07 mauro3

Unfortunately, it does not. I changed it and still had to kill the process after a minute. Doesn't look like the method is even going into that if clause

galbarel avatar Jul 01 '16 14:07 galbarel

This equation is stiff. Try a stiff solver (e.g. ode23s or DASSL.jl) instead, both can solve your problem very quickly.

EDIT: as a plan for a distant future we should add a stiffness detection system (sometimes it's difficult to see if the equation is stiff) and maybe switch between stiff and non stiff solvers.

pwl avatar Jul 01 '16 15:07 pwl

Also, equation by @galbarel is unrelated to the original issue (there is no backward time integration).

For a temporary remedy to the original issue you might want to take a look at #49 where backward time integration works (at least in the case @crbinz described).

pwl avatar Jul 01 '16 16:07 pwl

Thanks @pwl for noticing the stiff problem. The other problem seems to be something else, maybe related to backwards integration but not that as, note that e.g. ode78 or ODE.ode23 or ODE.ode45_fe work fine.

Also x-ref to the original julia-users thread https://groups.google.com/d/msg/julia-users/zB2mhbX3KeU/Tqn_hxWACQAJ .

mauro3 avatar Jul 01 '16 17:07 mauro3