Hello julia community,
I’m relatively new to julia and have a few questions regarding a dynamical system
representing a simple ecological model of five coupled equations.
The goal is to let the system integrate for some number of steps
and then look at qualities such as individual variances, wether it has been terminated
or the dominant lyapunov exponent of the last elements, repetitive for different parameters.
Here are my first steps towards the main aims:
using DifferentialEquations
using DynamicalSystems
using Statistics
function build(p1,p2)
function evolution!(du,u,p,t)
du[1] = u[1]*( 2.5*(1-u[1]) - (2.5*u[2])/(1+5*(u[1]+u[2]+u[3]) ) )
du[2] = u[2]*( 2.5*(1-u[2]) - (5*u[2])/(1+5*(u[1]+u[2]+u[3]) ) )
du[3] = u[3]*( 2.5*(1-u[3]) - (7.5*u[2])/(1+5*(u[1]+u[2]+u[3]) ) )
du[4] = u[4]*( ((2.5*u[1]+5*u[2]+7.5*u[3])/(1+5*(u[1]+u[2]+u[3]) )) - ((1*u[5])/(1+2*u[4])) - p[1] )
du[5] = u[5]*( (1*u[4])/(1+2*u[4]) - p[2] )
return nothing
end
function condition(u,t,integrator)
minimum(u)-10^-4
end
affect!(integrator) = terminate!(integrator)
cb = ContinuousCallback(condition,affect!)
u0 = SVector{5}(0.5,0.5,0.5,0.5,0.5)
ds = ContinuousDynamicalSystem(evolution!, u0, [0.1,0.1],t0 = 0.0)
for i in p1, j in p2
set_parameter!(ds,[i,j])
prob = ODEProblem(ds, (0,100.0),[i,j])
sol = solve(prob,callback=cb,dense=false)
if sol.t[end] != 100.0
println(sol.t[end])
elseif var(sol[1,end-20:end]) < 10^-4 && var(sol[2,end-20:end]) < 10^-4
println("equilibrium")
end
println(lyapunov(ds,100))
end
end
Calling
build([0.8],[0.4])
only takes a couple of seconds, but allready 8*8 different values seems never to end.
Now there a quite a few issues I got, any help to any of them is very much appreciated, even some link if I did’nt research enough. The integration period is aimed at a lot of steps of which the most intensive beside the integration I believe is going to be the computation of the lyapunov exponent of about the last fifth. That the whole process is going to take a while to compute is okay, I just want to optimize as much as possible.
i) I had the feeling some inplace static equations had roughly the same consumption for one solve, but is this now the right choice for so many iterations?
ii) I had the feeling that the trajectory() out of DynamicalSystems.jl had a slightly different behaviour, which I’d like to test but did’nt manage to implement the terminal callback with something like
trajectory(ds,100.0, callback = cb)
iii) Im relatively new to numerical integration and am rather insecure about the choice of solving algorithm, is the standard Tsit5() a good choice here? Is it reasonable to try out for higher tolerances explicitely? The numbers usually stay pretty low and I believe there are no huge slopes to expect.
iv)Do I have to specify some kind of multithreading or does julia automate the computation on the different cores? Would it be good to try out some form such that the solver is done one after another using the same space and storing the evaluation somewhere else because intermediate results are not required later? Regarding that which variables are better static like the initial conditions here for instance?
v) Is there property to check whether solver got terminated maybe by the callback, easier then comparing the final time-value?
vi) Is the usage of the minimum() a good choice in the callback to define some lower threshold or would it increase the performance by putting the callback aside and looking later through the values? The point is that in some cases after the threshold is reached another equation goes through the roof ~e80.
vii) If I wanted to implement more detailed conditional computation, say check for termination and only if not for the variance of all equations, are there faster ways then simple elseif?
viii) Is the
set_parameter!()
a bit expendable here as I already call the parameters in the ODEproblem?
Many Thanks for your time and effort people!