From 2c86306eca02251e1fe6b43d45833762264d9080 Mon Sep 17 00:00:00 2001 From: Torkel Date: Sat, 14 Oct 2023 11:18:38 -0400 Subject: [PATCH] updates --- Project.toml | 3 +- docs/Project.toml | 1 + docs/pages.jl | 1 + .../bifurcation_diagram_computation.md | 93 +++++++++++++++++++ ext/MTKBifurcationKitExt.jl | 18 ++-- test/extensions/bifurcationkit.jl | 29 ++++++ test/runtests.jl | 1 + 7 files changed, 136 insertions(+), 10 deletions(-) create mode 100644 docs/src/tutorials/bifurcation_diagram_computation.md create mode 100644 test/extensions/bifurcationkit.jl diff --git a/Project.toml b/Project.toml index ae899f1220..86da45a0a7 100644 --- a/Project.toml +++ b/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"] diff --git a/docs/Project.toml b/docs/Project.toml index 81a99e004e..68d33cc4c5 100644 --- a/docs/Project.toml +++ b/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" diff --git a/docs/pages.jl b/docs/pages.jl index 0b6d0555cd..0104c5eb48 100644 --- a/docs/pages.jl +++ b/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", 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 index 60ee91b9df..9a966ba55c 100644 --- a/ext/MTKBifurcationKitExt.jl +++ b/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 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 e6387e2cd9..0322bba523 100644 --- a/test/runtests.jl +++ b/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")