Merge pull request #2321 from TorkelE/BifurcationKitExt2

BifurcationKit Extension (take 2)

Project.toml

@@ -51,9 +51,11 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [weakdeps]
+BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6"
 
 [extensions]
+MTKBifurcationKitExt = "BifurcationKit"
 MTKDeepDiffsExt = "DeepDiffs"
 
 [compat]
@@ -101,6 +103,7 @@ julia = "1.6"
 [extras]
 AmplNLWriter = "7c4d4715-977e-5154-bfe0-e096adeac482"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 ControlSystemsMTK = "687d7614-c7e5-45fc-bfc3-9ee385575c88"
 DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
@@ -123,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", "BenchmarkTools", "BifurcationKit", "ControlSystemsMTK", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq"]

docs/pages.jl

@@ -8,6 +8,7 @@ pages = [
         "tutorials/programmatically_generating.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

@@ -0,0 +1,38 @@
+module MTKBifurcationKitExt
+
+### Preparations ###
+
+# Imports
+using ModelingToolkit, Setfield
+import BifurcationKit
+
+### Creates BifurcationProblem Overloads ###
+
+# When input is a NonlinearSystem.
+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=jac)
+    F = ofun.f
+    J = jac ? ofun.jac : nothing
+
+    # Computes bifurcation parameter and plot var indexes.
+    bif_idx = findfirst(isequal(bif_par), parameters(nsys))
+    if !isnothing(plot_var) 
+        plot_idx = findfirst(isequal(plot_var), states(nsys))
+        record_from_solution = (x, p) -> x[plot_idx]
+    end
+
+    # Converts the input state guess.
+    u0_bif = ModelingToolkit.varmap_to_vars(u0_bif, states(nsys))
+    ps = ModelingToolkit.varmap_to_vars(ps, parameters(nsys))
+
+    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, 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

@@ -59,3 +59,4 @@ end
 if VERSION >= v"1.9"
     @safetestset "Latexify recipes Test" include("latexify.jl")
 end
+@safetestset "BifurcationKit Extension Test" include("extensions/bifurcationkit.jl.jl")