How to catch integrator problems

Specific Domains Modelling & Simulations
1

I am using a CoupledODEs as a test for the results of a genetic algorithm. In other words, the parameter values (here 0.5 and 1.8) are chosen at random by the GA, then tested by solving the ODEs. This means that occasionally the parameter values chosen are a bit silly and lead to explosive results. In many cases I can catch these using various tests, but this one appears to go into an infinite loop, or at least, it doesn’t get caught by the try-block. Is there any way of catching and handling this?

using DynamicalSystems

function demo()
	odes = CoupledODEs(
		function (du,u,p,t)
			q = sum(p[1:2].*u[1:2])

			du[1] = (q/(q+u[3]) * (1 - p[1]/(p[1]+u[3])) - 0.1) * u[1]
			du[2] = (q/(q+u[3]) * (1 - p[2]/(p[2]+u[3])) - 0.1) * u[2]
			du[3] = 1 - u[3] - q
			nothing
		end,
		[2.0,2.0,1.0],
		[0.5,1.8]
	)

	for t in 0:0.1:0.4
        while t > current_time(odes)
			try
				step!(odes)
			catch (e)
				println("Caught!")
			end
        end
	end
end
2

Recently I was in exactly the same situation as you and CHris recommended me to add a dtmin option into the the ODE solver keywords. Something very small compared to the timescales of the system, e.g., 1e-3 or 1e-6 the typical timescale. What could happen is that the solver makes smaller and smaller steps.

Can you try this and let me know?

3

Notice that this try block doesn’t really help because the whole purpose of step! is to never explicitly throw an error, which is what the try blocks try to catch.

4

Thanks. Can you tell me where I would insert this dtmin option? I’m having some difficulty finding my way around the DynamicalSystems documentation - in particular regarding the step! function. When I click on the link to step! in the tutorial documentation, it generates a 404 File Not Found error. The code in ArbitrarySteppable gives me no clue, and indeed, in the above code, VSE underlines my call to step! as being a possible method call error. I get the feeling I’m doing something wrong somewhere. :thinking:

5

OK, sorry - think I’ve found what you mean in trajectory.jl. I’ll try it out.

6

Nope - I can’t work it out. I inserted the extra argument diffeq = (dtmin = 1e-4) after the parameters in the CoupledODEs constructor, but then I get the error:

MethodError: no method matching haskey(::Float64, ::Symbol)

generated from continuous_time_ode.jl:16

7

OK, finally got the diffeq issue sorted - finally realised I also had to specify the solver algorithm. Inserted this argument in the CoupledODEs() constructor:

diffeq = (alg = Tsit5(), dtmin = 1e-6)

But I’m afraid it still hangs at the same crucial time. However, the messages may be of use to you, so here they are:

┌ Warning: dt(1.0e-6) <= dtmin(1.0e-6) at t=0.38987336473510126, and step error estimate = 54.93179322897898. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase C:\Users\hswt136nia\.julia\packages\SciMLBase\eK30d\src\integrator_interface.jl:599