From 78bcb4eeb38334e1aa08c0c080c74d7984909dab Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 17 Oct 2023 21:26:30 +0100 Subject: [PATCH 1/3] init --- Project.toml | 5 +- docs/pages.jl | 1 + .../bifurcation_diagram_computation.md | 93 +++++++++++++++++++ ext/MTKBifurcationKitExt.jl | 40 ++++++++ test/extensions/bifurcationkit.jl | 29 ++++++ test/runtests.jl | 1 + 6 files changed, 168 insertions(+), 1 deletion(-) create mode 100644 docs/src/tutorials/bifurcation_diagram_computation.md create mode 100644 ext/MTKBifurcationKitExt.jl create mode 100644 test/extensions/bifurcationkit.jl diff --git a/Project.toml b/Project.toml index 4474aa35ad..446bfbeefe 100644 --- a/Project.toml +++ b/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"] diff --git a/docs/pages.jl b/docs/pages.jl index f6f5d7af50..9d49c477ec 100644 --- a/docs/pages.jl +++ b/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", diff --git a/docs/src/tutorials/bifurcation_diagram_computation.md b/docs/src/tutorials/bifurcation_diagram_computation.md new file mode 100644 index 0000000000..1bb6e33935 --- /dev/null +++ b/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. \ No newline at end of file diff --git a/ext/MTKBifurcationKitExt.jl b/ext/MTKBifurcationKitExt.jl new file mode 100644 index 0000000000..9a966ba55c --- /dev/null +++ b/ext/MTKBifurcationKitExt.jl @@ -0,0 +1,40 @@ +module MTKBifurcationKitExt + +println("BifurcationKit extension loaded") + +### 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 diff --git a/test/extensions/bifurcationkit.jl b/test/extensions/bifurcationkit.jl new file mode 100644 index 0000000000..c9ae8890ad --- /dev/null +++ b/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 \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index b60d82a6ce..10c56f3676 100644 --- a/test/runtests.jl +++ b/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") \ No newline at end of file From 65f41de87022a584f4cc8fab314fa18662922b8d Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 17 Oct 2023 21:33:31 +0100 Subject: [PATCH 2/3] fix --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 446bfbeefe..7bac7a4843 100644 --- a/Project.toml +++ b/Project.toml @@ -126,4 +126,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["AmplNLWriter", "BenchmarkTools", "BifurcationKit",, "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"] From eeb160b61a53f97012778ec94e782d2984c4ca10 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 17 Oct 2023 21:34:06 +0100 Subject: [PATCH 3/3] update --- ext/MTKBifurcationKitExt.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/ext/MTKBifurcationKitExt.jl b/ext/MTKBifurcationKitExt.jl index 9a966ba55c..2eb7b8807b 100644 --- a/ext/MTKBifurcationKitExt.jl +++ b/ext/MTKBifurcationKitExt.jl @@ -1,7 +1,5 @@ module MTKBifurcationKitExt -println("BifurcationKit extension loaded") - ### Preparations ### # Imports