-
-
Notifications
You must be signed in to change notification settings - Fork 210
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Bifurcation kit extension #2300
Closed
Closed
Changes from 3 commits
Commits
Show all changes
63 commits
Select commit
Hold shift + click to select a range
e56de20
split parameters linearize bug
cee69ed
added old MatrixGain definition for comparison
b4f821e
new promote_to_concrete
e62e3ee
improvements to protmote_to_concrete
ee4deb9
remake issue
a9c56b6
Suppose heterogeneous parameters for linearize and remake
YingboMa 1a96915
Fix zero_dummy_der
YingboMa 023e024
Relax tests
YingboMa 119cb27
Don't use SciMLBase 2
YingboMa 6bb0b2e
Merge branch 'master' into bgc/split_params_bug
YingboMa 29aec7b
Fix eq ordering
YingboMa addf335
Change bounds back
YingboMa 2e4c1bb
Format
YingboMa 2a553a7
Be careful with string sorting
YingboMa f9ea3a6
Update LaTeXify ref tests
YingboMa c2f92f1
init
TorkelE b05bc5f
makefunction
TorkelE f2d80b0
1.6 has a different term ordering
YingboMa ae05041
Fix DiscreteProblem construction
YingboMa 2309f9f
Merge pull request #2283 from SciML/bgc/split_params_bug
YingboMa fb7c3af
Update Project.toml
YingboMa 6e795b8
Fix recursive structure and setup precompilation
ChrisRackauckas 7949761
Merge pull request #2301 from SciML/precompile_workload
ChrisRackauckas c5bd821
docs(refactor): update acausal components tutorial with `@mtkmodel`
ven-k 6ac545d
docs(refactor): update structural identifiability tutorial with `@mtk…
ven-k 278698b
docs(fix): set a `u0` within the constraint for the optimization problem
ven-k 94fb749
docs(refactor): refactor "Getting Started" section with `@mtkmodel`
ven-k 438a225
docs(refactor): rename the `@mtkmodel` docs page as "Defining compone…
ven-k ff2a778
docs: reword Components and Connectors and elaborate Connectors and M…
ven-k 2ef4752
docs: add `connect` section to Variable metadata
ven-k 3034c9f
Merge pull request #2286 from ven-k/vkb/mtkmodel-examples
ChrisRackauckas c4868e3
Add `@mtkbuild` greatly simplify the tutorials
ChrisRackauckas ac75e77
Add `parent` to hierarchical systems and set it in `structural_simplify`
YingboMa b0c2938
Format
YingboMa c3c030e
Add "docs/src/assets/*.toml" to gitignore
YingboMa cd0f04a
Add tests
YingboMa 0d5ab59
Fix CI
YingboMa 7370a45
fix tutorial syntax
ChrisRackauckas 137d3d1
Support implicit name unpack in `at extend`
YingboMa c679dd9
fix acausal tutorial
ChrisRackauckas 17a3b6a
Merge pull request #2305 from SciML/mtkbuild
ChrisRackauckas 3d83c4c
one more tutorial simplification
ChrisRackauckas 0a9dca1
Merge branch 'master' into myb/extend
YingboMa 73f9e0a
Stop using `at testset` macro
YingboMa 718ae81
Remove ODAEProblem from acausal tutorial
ChrisRackauckas 0bcdee8
Merge pull request #2308 from SciML/odaeproblem
YingboMa 379161f
Merge pull request #2307 from SciML/tuts
YingboMa 314bb1d
Merge pull request #2306 from SciML/myb/extend
YingboMa 298a9a5
patch system names
ChrisRackauckas bdc45d1
Simplify `at extend`
YingboMa 1152f30
update docs for new extension syntax
ChrisRackauckas 9e45fee
fix for simplified extend
ChrisRackauckas 563f04f
Merge pull request #2310 from SciML/myb/extend
ChrisRackauckas 2c17d63
Merge pull request #2309 from SciML/extend
ChrisRackauckas d4239a7
Update Project.toml
ChrisRackauckas bc93796
non canonical docstrings
ArnoStrouwen b54b05b
update contributing docs
ArnoStrouwen 0ef0bb4
Merge pull request #2315 from ArnoStrouwen/docs
ChrisRackauckas 2c86306
updates
TorkelE 3249aeb
init
TorkelE 7f5d8be
makefunction
TorkelE c4a1a20
updates
TorkelE 918f5e6
Merge remote-tracking branch 'TorkelE/BifurcationKitExtension' into B…
TorkelE File filter
Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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. |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
fix the naming.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
what should I name it?