updates

Project.toml

@@ -102,6 +102,7 @@ julia = "1.6"
 
 [extras]
 AmplNLWriter = "7c4d4715-977e-5154-bfe0-e096adeac482"
+BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ControlSystemsMTK = "687d7614-c7e5-45fc-bfc3-9ee385575c88"
 DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6"
@@ -125,4 +126,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["AmplNLWriter", "BenchmarkTools", "ControlSystemsMTK", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq"]
+test = ["AmplNLWriter", "BifurcationKit", "BenchmarkTools", "ControlSystemsMTK", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq"]

docs/Project.toml

@@ -1,5 +1,6 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"

docs/pages.jl

@@ -7,6 +7,7 @@ pages = [
         "tutorials/modelingtoolkitize.md",
         "tutorials/stochastic_diffeq.md",
         "tutorials/parameter_identifiability.md",
+        "tutorials/bifurcation_diagram_computation.md",
         "tutorials/domain_connections.md"],
     "Examples" => Any["Basic Examples" => Any["examples/higher_order.md",
             "examples/spring_mass.md",

docs/src/tutorials/bifurcation_diagram_computation.md

@@ -0,0 +1,93 @@
+# [Bifurcation Diagrams](@id 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](https://github.com/bifurcationkit/BifurcationKit.jl) package. ModelingToolkit provides a simple interface for creating BifurcationKit compatible `BifurcationProblem`s from `NonlinearSystem`s and `ODESystem`s. All teh 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](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/).
+
+### Creating a `BifurcationProblem`
+Let us first consider a simple `NonlinearSystem`:
+```@example Bif1
+using ModelingToolkit
+@variables t x(t) y(t)
+@parameters μ α
+eqs = [0 ~ μ*x - x^3 + α*y,
+        0 ~ -y]
+@named nsys = NonlinearSystem(eqs, [x, y], [μ, α])
+```
+we wish to compute a bifurcation diagram for this system as we vary the parameter `μ`. For this, we need to provide the following information:
+1. The system for which we wish to compute the bifurcation diagram (`nsys`).
+2. The parameter which we wish to vary (`μ`).
+3. The parameter set for which we want to compute the bifurcation diagram.
+4. An initial guess of the state of the system for which there is a steady state at our provided parameter value.
+5. 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:
+```@example Bif1
+bif_par = μ
+p_start = [μ => -1.0, α => 1.0]
+u0_guess = [x => 1.0, y => 1.0]
+plot_var = x;
+```
+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.
+```@example Bif1
+using BifurcationKit
+bprob = BifurcationProblem(nsys, u0_guess, p_start, bif_par; plot_var=plot_var, jac=false)
+```
+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). 
+```@example Bif1
+p_span = (-4.0, 6.0)
+opt_newton = NewtonPar(tol = 1e-9, max_iterations = 20)
+opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds = 0.01,
+	max_steps = 100, nev = 2, newton_options = opt_newton,
+	p_min = p_span[1], p_max = p_span[2],
+	detect_bifurcation = 3, n_inversion = 4, tol_bisection_eigenvalue = 1e-8, dsmin_bisection = 1e-9);
+```
+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.
+```@example Bif1
+bf = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside=true)
+```
+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.
+```@example Bif1
+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 `ODESystem`s (rather than `NonlinearSystem`s) as input to `BifurcationProblem`. Here follows a brief such example.
+
+```@example Bif2
+using BifurcationKit, ModelingToolkit, Plots
+
+@variables t x(t) y(t)
+@parameters μ
+D = Differential(t)
+eqs = [D(x) ~ μ*x - y - x*(x^2+y^2), 
+        D(y) ~ x + μ*y - y*(x^2+y^2)]
+@named 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)
+opt_newton = NewtonPar(tol = 1e-9, max_iterations = 20)
+opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds = 0.01,
+	max_steps = 100, nev = 2, newton_options = opt_newton,
+	p_max = p_span[2], p_min = p_span[1],
+	detect_bifurcation = 3, n_inversion = 4, tol_bisection_eigenvalue = 1e-8, dsmin_bisection = 1e-9)
+
+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.

ext/MTKBifurcationKitExt.jl

@@ -11,11 +11,11 @@ import BifurcationKit
 ### Creates BifurcationProblem Overloads ###
 
 # When input is a NonlinearSystem.
-function BifurcationKit.BifurcationProblem(nsys::NonlinearSystem, u0_guess, p_start, bif_par, args...; plot_var=nothing, record_from_solution=BifurcationKit.record_sol_default, kwargs...)
+function BifurcationKit.BifurcationProblem(nsys::NonlinearSystem, u0_bif, ps, bif_par, args...; plot_var=nothing, record_from_solution=BifurcationKit.record_sol_default, jac=true, kwargs...)
     # Creates F and J functions.
-    ofun = NonlinearFunction(nsys; jac=true)
+    ofun = NonlinearFunction(nsys; jac=jac)
     F = ofun.f
-    J = ofun.jac
+    J = jac ? ofun.jac : nothing
 
     # Computes bifurcation parameter and plot var indexes.
     bif_idx = findfirst(isequal(bif_par), parameters(nsys))
@@ -25,16 +25,16 @@ function BifurcationKit.BifurcationProblem(nsys::NonlinearSystem, u0_guess, p_st
     end
 
     # Converts the input state guess.
-    u0_guess = ModelingToolkit.varmap_to_vars(u0_guess, states(nsys))
+    u0_bif = ModelingToolkit.varmap_to_vars(u0_bif, states(nsys))
+    ps = ModelingToolkit.varmap_to_vars(ps, parameters(nsys))
 
-    return BifurcationKit.BifurcationProblem(F, u0_guess, [p_start], (@lens _[1]), args...; record_from_solution = record_from_solution, J = J, kwargs...)
+    return BifurcationKit.BifurcationProblem(F, u0_bif, ps, (@lens _[bif_idx]), args...; record_from_solution = record_from_solution, J = J, kwargs...)
 end
 
 # When input is a ODESystem.
-function BifurcationKit.BifurcationProblem(osys::ODESystem, u0, p, args...; kwargs...)
-    return BifurcationProblem(convert(NonlinearProblem, osys), args...; kwargs...)
+function BifurcationKit.BifurcationProblem(osys::ODESystem, args...; kwargs...)
+    nsys = NonlinearSystem([0 ~ eq.rhs for eq in equations(osys)], states(osys), parameters(osys); name=osys.name)
+    return BifurcationKit.BifurcationProblem(nsys, args...; kwargs...)
 end
 
-
-
 end # module

test/extensions/bifurcationkit.jl

@@ -0,0 +1,29 @@
+using BifurcationKit, ModelingToolkit, Test
+
+# Checks pitchfork diagram and that there are the correct number of branches (a main one and two children)
+let 
+    @variables t x(t) y(t)
+    @parameters μ α
+    eqs = [0 ~ μ*x - x^3 + α*y,
+            0 ~ -y]
+    @named nsys = NonlinearSystem(eqs, [x, y], [μ, α])
+
+    bif_par = μ
+    p_start = [μ => -1.0, α => 1.0]
+    u0_guess = [x => 1.0, y => 1.0]
+    plot_var = x;
+
+    using BifurcationKit
+    bprob = BifurcationProblem(nsys, u0_guess, p_start, bif_par; plot_var=plot_var, jac=false)
+
+    p_span = (-4.0, 6.0)
+    opt_newton = NewtonPar(tol = 1e-9, max_iterations = 20)
+    opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds = 0.01,
+        max_steps = 100, nev = 2, newton_options = opt_newton,
+        p_min = p_span[1], p_max = p_span[2],
+        detect_bifurcation = 3, n_inversion = 4, tol_bisection_eigenvalue = 1e-8, dsmin_bisection = 1e-9);
+
+    bf = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside=true)
+
+    @test length(bf.child) == 2
+end

test/runtests.jl

@@ -55,5 +55,6 @@ end
 @safetestset "OptimizationSystem Test" include("optimizationsystem.jl")
 @safetestset "FuncAffect Test" include("funcaffect.jl")
 @safetestset "Constants Test" include("constants.jl")
+@safetestset "BifurcationKit Extension Test" include("extensions/bifurcationkit.jl")
 # Reference tests go Last
 @safetestset "Latexify recipes Test" include("latexify.jl")