Skip to content

Commit

Permalink
fix: BifurcationKit test path + format
Browse files Browse the repository at this point in the history
  • Loading branch information
ven-k committed Oct 19, 2023
1 parent 31ed28e commit 442cba7
Show file tree
Hide file tree
Showing 7 changed files with 104 additions and 39 deletions.
1 change: 1 addition & 0 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
1 change: 1 addition & 0 deletions docs/src/systems/JumpSystem.md
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ structural_simplify
```@docs; canonical=false
DiscreteProblem(sys::ModelingToolkit.DiscreteSystem, u0map, tspan)
```

```@docs
JumpProblem(sys::JumpSystem, prob, aggregator)
```
1 change: 1 addition & 0 deletions docs/src/systems/SDESystem.md
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ sde = SDESystem(ode, noiseeqs)
structural_simplify
alias_elimination
```

```@docs
ModelingToolkit.Girsanov_transform
```
Expand Down
88 changes: 63 additions & 25 deletions docs/src/tutorials/bifurcation_diagram_computation.md
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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.

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.
28 changes: 23 additions & 5 deletions ext/MTKBifurcationKitExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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,

Check warning on line 12 in ext/MTKBifurcationKitExt.jl

View check run for this annotation

Codecov / codecov/patch

ext/MTKBifurcationKitExt.jl#L12

Added line #L12 was not covered by tests
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)

Check warning on line 22 in ext/MTKBifurcationKitExt.jl

View check run for this annotation

Codecov / codecov/patch

ext/MTKBifurcationKitExt.jl#L22

Added line #L22 was not covered by tests
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)

Check warning on line 28 in ext/MTKBifurcationKitExt.jl

View check run for this annotation

Codecov / codecov/patch

ext/MTKBifurcationKitExt.jl#L28

Added line #L28 was not covered by tests
plot_idx = findfirst(isequal(plot_var), states(nsys))
record_from_solution = (x, p) -> x[plot_idx]
end
Expand All @@ -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,

Check warning on line 37 in ext/MTKBifurcationKitExt.jl

View check run for this annotation

Codecov / codecov/patch

ext/MTKBifurcationKitExt.jl#L37

Added line #L37 was not covered by tests
u0_bif,
ps,
(@lens _[bif_idx]),

Check warning on line 40 in ext/MTKBifurcationKitExt.jl

View check run for this annotation

Codecov / codecov/patch

ext/MTKBifurcationKitExt.jl#L40

Added line #L40 was not covered by tests
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)],

Check warning on line 49 in ext/MTKBifurcationKitExt.jl

View check run for this annotation

Codecov / codecov/patch

ext/MTKBifurcationKitExt.jl#L49

Added line #L49 was not covered by tests
states(osys),
parameters(osys);
name = osys.name)
return BifurcationKit.BifurcationProblem(nsys, args...; kwargs...)
end

Expand Down
22 changes: 14 additions & 8 deletions test/extensions/bifurcationkit.jl
Original file line number Diff line number Diff line change
@@ -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
end
2 changes: 1 addition & 1 deletion test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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")
@safetestset "BifurcationKit Extension Test" include("extensions/bifurcationkit.jl")

0 comments on commit 442cba7

Please sign in to comment.