# Plot model stability region depending on parameters

**URL:** <https://discourse.julialang.org/t/plot-model-stability-region-depending-on-parameters/79584>\
**Category:** Modelling & Simulations\
**Tags:** question\
**Created:** [April 16, 2022, 11:51pm UTC](https://discourse.julialang.org/t/plot-model-stability-region-depending-on-parameters/79584 "2022-04-16T23:51:52Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![Sciemon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sciemon/32/20897_2.png) [@Sciemon](https://discourse.julialang.org/u/Sciemon)\
**Post date:** [April 16, 2022, 11:51pm UTC](https://discourse.julialang.org/t/plot-model-stability-region-depending-on-parameters/79584/1 "2022-04-16T23:51:53Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 17, 2022, 6:29am UTC](https://discourse.julialang.org/t/plot-model-stability-region-depending-on-parameters/79584/2 "2022-04-17T06:29:10Z")

</div>

What do you mean by `A` and `B` here? Are both part of the Jacobian?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 17, 2022, 6:33am UTC](https://discourse.julialang.org/t/plot-model-stability-region-depending-on-parameters/79584/3 "2022-04-17T06:33:15Z")

</div>

You could likely make use of  
[https://docs.juliaplots.org/stable/gallery/gr/generated/gr-ref22/#gr\_demo\_22](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.

---

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [April 17, 2022, 7:36am UTC](https://discourse.julialang.org/t/plot-model-stability-region-depending-on-parameters/79584/4 "2022-04-17T07:36:49Z")

</div>

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:

```julia
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)`.

---

<div class="post-metadata">

**Author:** ![Sciemon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sciemon/32/20897_2.png) [@Sciemon](https://discourse.julialang.org/u/Sciemon)\
**Post date:** [April 17, 2022, 4:25pm UTC](https://discourse.julialang.org/t/plot-model-stability-region-depending-on-parameters/79584/5 "2022-04-17T16:25:03Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Sciemon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sciemon/32/20897_2.png) [@Sciemon](https://discourse.julialang.org/u/Sciemon)\
**Post date:** [April 17, 2022, 4:34pm UTC](https://discourse.julialang.org/t/plot-model-stability-region-depending-on-parameters/79584/6 "2022-04-17T16:34:10Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [May 1, 2022, 11:22am UTC](https://discourse.julialang.org/t/plot-model-stability-region-depending-on-parameters/79584/7 "2022-05-01T11:22:02Z")

</div>

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:

- [fold curve](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/ode/tutorialCO/#CO-oxydation-(codim-2))
- [richer example](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/ode/lorenz84/#Extended-Lorenz-84-model-(codim-2-BT/ZH-aBS))

The beauty of this is that you dont need to sample the 2d space (A,B) 😉

---

<div class="post-metadata">

**Author:** ![Sciemon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sciemon/32/20897_2.png) [@Sciemon](https://discourse.julialang.org/u/Sciemon)\
**Post date:** [May 3, 2022, 8:56am UTC](https://discourse.julialang.org/t/plot-model-stability-region-depending-on-parameters/79584/8 "2022-05-03T08:56:01Z")

</div>

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”?

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [May 3, 2022, 10:23am UTC](https://discourse.julialang.org/t/plot-model-stability-region-depending-on-parameters/79584/9 "2022-05-03T10:23:27Z")

</div>

Kuznetsov, Y. A. _Elements of Applied Bifurcation Theory_ . New York, NY: Springer New York, 2004. [Elements of Applied Bifurcation Theory | SpringerLink](http://dx.doi.org/10.1007/978-1-4757-3978-7).

I guess I could also improve [docs](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/guidelines/) to better explain this