Plot model stability region depending on parameters

Specific Domains Modelling & Simulations
4

I’m certain you can also do this with DiffEq.jl, because they also calculate a jacobian, but here is the DynamicalSystems.jl version that I find handy:

using DynamicalSystems, GLMakie, LinearAlgebra
isstable(e) = maximum(real(x) for x in e) < 0
function dynamic_rule(u, p, t)
A, B = p
# insert your actual equations here 
# (not linearized)
end

u0 = [ ] # point you are interested in
ds = ContinuousDynamicalSystem(dynamic_rule, u0, p)

stab_matrix = zeros(length(As), length(Bs))
for (i, A) in enumerate(As)
set_parameter!(ds, 1, A)
for (j, B) in enumerate(Bs)
set_parameter!(ds, 2, B)
J = jacobian(ds)
e = LinearAlgebra.eigvals(Array(J))
stab_matrix[i,j] = isstable(e)
end
end

heatmap(stab_matrix)

BTW if you already have an ODEProblem you can just do ds = ContinuousDynamicalSystem(odeprob).

1 Like