From 442cba7613ca42872654cfef404e5ec5fb624fa9 Mon Sep 17 00:00:00 2001 From: Venkateshprasad <32921645+ven-k@users.noreply.github.com> Date: Thu, 19 Oct 2023 15:40:00 +0530 Subject: [PATCH] fix: BifurcationKit test path + format --- Project.toml | 1 + docs/src/systems/JumpSystem.md | 1 + docs/src/systems/SDESystem.md | 1 + .../bifurcation_diagram_computation.md | 88 +++++++++++++------ ext/MTKBifurcationKitExt.jl | 28 ++++-- test/extensions/bifurcationkit.jl | 22 +++-- test/runtests.jl | 2 +- 7 files changed, 104 insertions(+), 39 deletions(-) diff --git a/Project.toml b/Project.toml index 7bac7a4843..9426118ed1 100644 --- a/Project.toml +++ b/Project.toml @@ -33,6 +33,7 @@ MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" +PrettyPrint = "8162dcfd-2161-5ef2-ae6c-7681170c5f98" RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47" diff --git a/docs/src/systems/JumpSystem.md b/docs/src/systems/JumpSystem.md index 68f0de2a5c..b066dbfc44 100644 --- a/docs/src/systems/JumpSystem.md +++ b/docs/src/systems/JumpSystem.md @@ -26,6 +26,7 @@ structural_simplify ```@docs; canonical=false DiscreteProblem(sys::ModelingToolkit.DiscreteSystem, u0map, tspan) ``` + ```@docs JumpProblem(sys::JumpSystem, prob, aggregator) ``` diff --git a/docs/src/systems/SDESystem.md b/docs/src/systems/SDESystem.md index 17f062c5ae..5e4dc958c0 100644 --- a/docs/src/systems/SDESystem.md +++ b/docs/src/systems/SDESystem.md @@ -26,6 +26,7 @@ sde = SDESystem(ode, noiseeqs) structural_simplify alias_elimination ``` + ```@docs ModelingToolkit.Girsanov_transform ``` diff --git a/docs/src/tutorials/bifurcation_diagram_computation.md b/docs/src/tutorials/bifurcation_diagram_computation.md index 1bb6e33935..45dc98190b 100644 --- a/docs/src/tutorials/bifurcation_diagram_computation.md +++ b/docs/src/tutorials/bifurcation_diagram_computation.md @@ -1,65 +1,91 @@ # [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] +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). + + 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) +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. +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). +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); + 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) +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") +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 @@ -68,26 +94,38 @@ 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)] +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] +u0_guess = [x => 0.0, y => 0.0] -bprob = BifurcationProblem(osys, u0_guess, p_start, bif_par; plot_var=plot_var, jac=false) +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) + 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) +bf = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside = true) using Plots -plot(bf; putspecialptlegend=false, markersize=2, plotfold=false, xguide="μ", yguide = "x") +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 + +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. diff --git a/ext/MTKBifurcationKitExt.jl b/ext/MTKBifurcationKitExt.jl index 2eb7b8807b..c85011b65f 100644 --- a/ext/MTKBifurcationKitExt.jl +++ b/ext/MTKBifurcationKitExt.jl @@ -9,15 +9,23 @@ 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...) +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) + 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) + if !isnothing(plot_var) plot_idx = findfirst(isequal(plot_var), states(nsys)) record_from_solution = (x, p) -> x[plot_idx] end @@ -26,12 +34,22 @@ function BifurcationKit.BifurcationProblem(nsys::NonlinearSystem, u0_bif, ps, bi 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...) + 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) + nsys = NonlinearSystem([0 ~ eq.rhs for eq in equations(osys)], + states(osys), + parameters(osys); + name = osys.name) return BifurcationKit.BifurcationProblem(nsys, args...; kwargs...) end diff --git a/test/extensions/bifurcationkit.jl b/test/extensions/bifurcationkit.jl index c9ae8890ad..45f0a89d57 100644 --- a/test/extensions/bifurcationkit.jl +++ b/test/extensions/bifurcationkit.jl @@ -1,29 +1,35 @@ using BifurcationKit, ModelingToolkit, Test # Checks pitchfork diagram and that there are the correct number of branches (a main one and two children) -let +let @variables t x(t) y(t) @parameters μ α - eqs = [0 ~ μ*x - x^3 + α*y, - 0 ~ -y] + 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; + plot_var = x using BifurcationKit - bprob = BifurcationProblem(nsys, u0_guess, p_start, bif_par; plot_var=plot_var, jac=false) + 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); + 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) + bf = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside = true) @test length(bf.child) == 2 -end \ No newline at end of file +end diff --git a/test/runtests.jl b/test/runtests.jl index 7d3ea5df68..5cafa89f68 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -59,4 +59,4 @@ end if VERSION >= v"1.9" @safetestset "Latexify recipes Test" include("latexify.jl") end -@safetestset "BifurcationKit Extension Test" include("extensions/bifurcationkit.jl") \ No newline at end of file +@safetestset "BifurcationKit Extension Test" include("extensions/bifurcationkit.jl")