Bifurcation Diagrams

Bifurcation diagrams describes how, for a dynamic system, the quantity and quality of its steady states changes with a parameter's value. These can be computed through the BifurcationKit.jl package. ModelingToolkit provides a simple interface for creating BifurcationKit compatible BifurcationProblems from NonlinearSystems and ODESystems. All the features provided by BifurcationKit can then be applied to these systems. This tutorial provides a brief introduction for these features, with BifurcationKit.jl providing a more extensive documentation.

Creating a `BifurcationProblem`

Let us first consider a simple NonlinearSystem:

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D

@variables x(t) y(t)
@parameters μ α
eqs = [0 ~ μ * x - x^3 + α * y,
    0 ~ -y]
@mtkbuild nsys = NonlinearSystem(eqs, [x, y], [μ, α])

\[ \begin{align} 0 &= x\left( t \right) \mu + y\left( t \right) \alpha - \left( x\left( t \right) \right)^{3} \end{align} \]

we wish to compute a bifurcation diagram for this system as we vary the parameter μ. For this, we need to provide the following information:

The system for which we wish to compute the bifurcation diagram (nsys).
The parameter which we wish to vary (μ).
The parameter set for which we want to compute the bifurcation diagram.
An initial guess of the state of the system for which there is a steady state at our provided parameter value.
The variable which value we wish to plot in the bifurcation diagram (this argument is optional, if not provided, BifurcationKit default plot functions are used).

We declare this additional information:

bif_par = μ
p_start = [μ => -1.0, α => 1.0]
u0_guess = [x => 1.0, y => 1.0]
plot_var = x;

\[ \begin{equation} x\left( t \right) \end{equation} \]

For the initial state guess (u0_guess), typically any value can be provided, however, read BifurcatioKit's documentation for more details.

We can now create our BifurcationProblem, which can be provided as input to BifurcationKit's various functions.

using BifurcationKit
bprob = BifurcationProblem(nsys,
    u0_guess,
    p_start,
    bif_par;
    plot_var = plot_var,
    jac = false)

┌─ Bifurcation Problem with uType Vector{Float64}
├─ Inplace:  true
├─ Dimension:  1
├─ Symmetric: false
└─ Parameter: p2

Here, the jac argument (by default set to true) sets whenever to provide BifurcationKit with a Jacobian or not.

Computing a bifurcation diagram

Let us consider the BifurcationProblem from the last section. If we wish to compute the corresponding bifurcation diagram we must first declare various settings used by BifurcationKit to compute the diagram. These are stored in a ContinuationPar structure (which also contain a NewtonPar structure).

p_span = (-4.0, 6.0)
opts_br = ContinuationPar(nev = 2,
    p_min = p_span[1],
    p_max = p_span[2])

BifurcationKit.ContinuationPar{Float64, BifurcationKit.DefaultLS, BifurcationKit.DefaultEig{typeof(real)}}
  dsmin: Float64 0.0001
  dsmax: Float64 0.1
  ds: Float64 0.01
  a: Float64 0.5
  p_min: Float64 -4.0
  p_max: Float64 6.0
  max_steps: Int64 400
  newton_options: BifurcationKit.NewtonPar{Float64, BifurcationKit.DefaultLS, BifurcationKit.DefaultEig{typeof(real)}}
  η: Float64 150.0
  save_to_file: Bool false
  save_sol_every_step: Int64 1
  nev: Int64 2
  save_eig_every_step: Int64 1
  save_eigenvectors: Bool true
  plot_every_step: Int64 10
  tol_stability: Float64 1.0e-10
  detect_fold: Bool true
  detect_bifurcation: Int64 3
  dsmin_bisection: Float64 1.0e-16
  n_inversion: Int64 2
  max_bisection_steps: Int64 25
  tol_bisection_eigenvalue: Float64 1.0e-16
  detect_event: Int64 0
  tol_param_bisection_event: Float64 1.0e-16
  detect_loop: Bool false

Here, p_span sets the interval over which we wish to compute the diagram.

Next, we can use this as input to our bifurcation diagram, and then plot it.

bf = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside = true)

[Bifurcation diagram]
 ┌─ From 0-th bifurcation point.
 ├─ Children number: 2
 └─ Root (recursion level 1)
      ┌─ Curve type: EquilibriumCont
      ├─ Number of points: 82
      ├─ Type of vectors: Vector{Float64}
      ├─ Parameter p2 starts at -4.0, ends at 6.0
      ├─ Algo: PALC
      └─ Special points:

- #  1, endpoint at p2 ≈ -4.00000000,                                                                     step =   0
- #  2,       bp at p2 ≈ +0.00044561 ∈ (-0.00010682, +0.00044561), |δp|=6e-04, [converged], δ = ( 1,  0), step =  38
- #  3, endpoint at p2 ≈ +6.00000000,                                                                     step =  81

Here, the value 2 sets how sub-branches of the diagram that BifurcationKit should compute. Generally, for bifurcation diagrams, it is recommended to use the bothside=true argument.

using Plots
plot(bf;
    putspecialptlegend = false,
    markersize = 2,
    plotfold = false,
    xguide = "μ",
    yguide = "x")

Here, the system exhibits a pitchfork bifurcation at μ=0.0.

Using `ODESystem` inputs

It is also possible to use ODESystems (rather than NonlinearSystems) as input to BifurcationProblem. Here follows a brief such example.

using BifurcationKit, ModelingToolkit, Plots
using ModelingToolkit: t_nounits as t, D_nounits as D

@variables x(t) y(t)
@parameters μ
eqs = [D(x) ~ μ * x - y - x * (x^2 + y^2),
    D(y) ~ x + μ * y - y * (x^2 + y^2)]
@mtkbuild osys = ODESystem(eqs, t)

bif_par = μ
plot_var = x
p_start = [μ => 1.0]
u0_guess = [x => 0.0, y => 0.0]

bprob = BifurcationProblem(osys,
    u0_guess,
    p_start,
    bif_par;
    plot_var = plot_var,
    jac = false)

p_span = (-3.0, 3.0)
opts_br = ContinuationPar(nev = 2,
    p_max = p_span[2], p_min = p_span[1])

bf = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside = true)
using Plots
plot(bf;
    putspecialptlegend = false,
    markersize = 2,
    plotfold = false,
    xguide = "μ",
    yguide = "x")

Here, the value of x in the steady state does not change, however, at μ=0 a Hopf bifurcation occur and the steady state turn unstable.

We compute the branch of periodic orbits which is nearby the Hopf Bifurcation. We thus provide the branch bf.γ, the index of the Hopf point we want to branch from: 2 in this case and a method PeriodicOrbitOCollProblem(20, 5) to compute periodic orbits.

br_po = continuation(bf.γ, 2, opts_br,
    PeriodicOrbitOCollProblem(20, 5);)

plot(bf; putspecialptlegend = false,
    markersize = 2,
    plotfold = false,
    xguide = "μ",
    yguide = "x")
plot!(br_po, xguide = "μ", yguide = "x", label = "Maximum of periodic orbit")

Let's see how to plot the periodic solution we just computed:

sol = get_periodic_orbit(br_po, 10)
plot(sol.t, sol[1, :], yguide = "x", xguide = "time", label = "")

Bifurcation Diagrams

Creating a BifurcationProblem

Computing a bifurcation diagram

Using ODESystem inputs

Creating a `BifurcationProblem`

Using `ODESystem` inputs