ode23s stalling

[jmaces]

When I run

ode23s(...)

on a linear ODE with ill-conditioned coefficient matrix, the estimated error err in line 302 returns NaN.

This results in being stuck in an infinite loop, since neither the current time t nor the step size h are changed in this case.

[acroy]

Would you mind to post an example? On Fr., 6. Nov. 2015 at 20:21, jmaces notifications@github.com wrote:

When I run

ode23s(...)

on a linear ODE with ill-conditioned coefficient matrix, the estimated error err in line 302 returns NaN.

This results in being stuck in an infinite loop, since neither the current time t nor the step size h are changed in this case.

— Reply to this email directly or view it on GitHub https://github.com/JuliaLang/ODE.jl/issues/81.

[jmaces]

Sure. I tried to simplify my original problem as much as possible and came up with an example that looks very harmless. However ode23s does not terminate.

using ODE

c = -120.0
k = 30.0

A = [0.0 1.0;-k -c]

#some info print outs
println(@sprintf("Conditioning: \t %f",cond(A)))
λ1 = 1/2*(-sqrt(c^2-4*k+0im)-c)
λ2 = 1/2*(sqrt(c^2-4*k+0im)-c)
Q = max(abs(real(λ1)),abs(real(λ2)))/min(abs(real(λ1)),abs(real(λ2)))
println(@sprintf("Eigenvalues: \t %f + %fi and %f + %fi",real(λ1),imag(λ1),real(λ2),imag(λ2)))
println(@sprintf("Stiffness: \t %f",Q))

#this line stalls the program
tout, yout = ode23s((t,y)->A*y,[10.0; 0.0], linspace(0.0,10.0,100))

I ran my tests on a JuliaBox server with an 0.3.11 kernel.

[acroy]

Thanks. I guess the easiest solution is checking isnan(err) and breaking out of the loop with an error message. Are you interested in submitting a PR with the fix & your test case? On Sa., 7. Nov. 2015 at 04:48, jmaces notifications@github.com wrote:

Sure. I tried to simplify my original problem as much as possible and came up with an example that looks very harmless. However ode32s does not terminate.

using ODE

c = -120.0 k = 30.0

A = [0.0 1.0;-k -c]

some info print outs

println(@sprintf("Conditioning: \t %f",cond(A))) λ1 = 1/2_(-sqrt(c^2-4k+0im)-c) λ2 = 1/2(sqrt(c^2-4_k+0im)-c) Q = max(abs(real(λ1)),abs(real(λ2)))/min(abs(real(λ1)),abs(real(λ2))) println(@sprintf("Eigenvalues: \t %f + %fi and %f + %fi",real(λ1),imag(λ1),real(λ2),imag(λ2))) println(@sprintf("Stiffness: \t %f",Q))

this line stalls the program

tout, yout = ode23s((t,y)->A*y,[10.0; 0.0], linspace(0.0,10.0,100))

I ran my tests on a JuliaBox server with an 0.3.11 kernel.

— Reply to this email directly or view it on GitHub https://github.com/JuliaLang/ODE.jl/issues/81#issuecomment-154614238.

[jmaces]

I can do that. Won't have time to do it until after the weekend though.

[mauro3]

The Runge Kutta solvers do this using the function isoutofdomain: https://github.com/JuliaLang/ODE.jl/blob/f0b83efacf7355b16adca96092fee33c17579443/src/ODE.jl#L91. Probably best if you also use that interface with ode23s. Also have a look at https://github.com/JuliaLang/ODE.jl/blob/9865c128fad90b9705c5aac5729bf456a93dc762/test/interface-tests.jl

[acroy]

@mauro3 : Good to know. We should probably add this to the API specs. I noticed that the RK solvers give a "td < minstep" error for @jmaces example. Wouldn't it be better to indicate that the solution went outside the domain?

[mauro3]

Yes, that would be better. I'll have a look at it. I should probably have a look at the API specs. I was reasonably thorough with documenting the API in the top comment of each function, but that should make it to the API specs.

[mauro3]

I had a look at it: outofdomain is just used inside the step-control. If it returns true, the step size is reduced, until it reaches the minimum. Then that error is thrown. I think, in general, it would be helpful to be more specific on what caused the minimum stepsize to be reached, but I don't have time to look into it at the moment.