I would like to plot the stability region of a (linearized) dynamic model at a certain operating point depending on model parameters. For example, assume that a model has parameters A and B. What I want is a plot that shows A on the x-axis and B on the y-axis and describes the region where the model is stable (or a line which describes the boundary between stability and instability). Instability would be defined as positive eigenvalue real parts of the Jacobian for a parameter pair. Is there a way to do this in Julia (starting from a model as defined for using DifferentialEquations.jl)? I looked at BifurcationKit.jl but couldn’t find an example that fits. Goes without saying I’m not the expert.
What do you mean by A and B here? Are both part of the Jacobian?
You could likely make use of
https://docs.juliaplots.org/stable/gallery/gr/generated/gr-ref22/#gr_demo_22
and apply it on a function that extracts the largest real value of the eigenvalues of the Jacobian. The values for the contours can be chosen manually, so you can make sure that zero is among (or the only) the contours.
1 Like
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
Thanks a lot for the immediate response! Yes this works and it is the brute force approach with many jacobian evaluations. I thought that maybe someone knows some magic how to do this more efficiently (maybe people that implement contour plots).
I once did this in Matlab and it was slow and I tried to reduce the number of jacobian evaluations by stepping along the boundary (imagine stepping between stable and unstable region in small steps to draw the boundary line). It was simple and worked quite well and drastically reduced the number of jacobian evaluations.
But for now I’m happy with this solution as Julia is generally faster.
Thanks, this is a great solution to draw the line! I used CairoMakie and implemented this as follows:
CairoMakie.contour(stab_matrix,levels = [0.0],colormap = [:red])
The stab_matrix contains the maximum eigenvalue real part at each point (looks smoother than using boolean values for stable/unstable). I define a colormap with only one color which sets the color of the stability boundary line.
Sounds like you are trying to plot codim 2 bifurcations. Say you detected a bifurcation for A fixed on the B axis. Can you trace this singularity in the A-B plane?
The answer is yes. Basic examples are:
The beauty of this is that you dont need to sample the 2d space (A,B) 
1 Like
Thanks! This looks very much like what I was looking for and will give it a try once I find time. As mentioned, I looked at BifurcationKit.jl, but couldn’t find the right bifurcation or right example for me. Are you aware of any resources (online and free at best 
) that provide an overview of the field for people who have never heard of “codim 2”?
Kuznetsov, Y. A. Elements of Applied Bifurcation Theory . New York, NY: Springer New York, 2004. Elements of Applied Bifurcation Theory | SpringerLink.
I guess I could also improve docs to better explain this
1 Like