Skip to content

Commit 50ff0fe

Browse files
Merge pull request #2445 from rveltz/patch-2
Update bifurcation_diagram_computation.md, add periodic orbits comput…
2 parents 3c756ee + 887f698 commit 50ff0fe

File tree

1 file changed

+25
-5
lines changed

1 file changed

+25
-5
lines changed

docs/src/tutorials/bifurcation_diagram_computation.md

Lines changed: 25 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -110,12 +110,9 @@ bprob = BifurcationProblem(osys,
110110
jac = false)
111111
112112
p_span = (-3.0, 3.0)
113-
opt_newton = NewtonPar(tol = 1e-9, max_iterations = 20)
114-
opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds = 0.01,
115-
max_steps = 100, nev = 2, newton_options = opt_newton,
113+
opts_br = ContinuationPar(nev = 2,
116114
p_max = p_span[2], p_min = p_span[1],
117-
detect_bifurcation = 3, n_inversion = 4, tol_bisection_eigenvalue = 1e-8,
118-
dsmin_bisection = 1e-9)
115+
)
119116
120117
bf = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside = true)
121118
using Plots
@@ -128,3 +125,26 @@ plot(bf;
128125
```
129126

130127
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.
128+
129+
We compute the branch of periodic orbits which is nearby the Hopf Bifurcation. We thus provide the branch `bf.γ`, the index of the Hopf point we want to branch from: 2 in this case and a method `PeriodicOrbitOCollProblem(20, 5)` to compute periodic orbits.
130+
131+
```@example Bif2
132+
br_po = continuation(bf.γ, 2, opts_br,
133+
PeriodicOrbitOCollProblem(20, 5);
134+
)
135+
136+
plot(bf; putspecialptlegend = false,
137+
markersize = 2,
138+
plotfold = false,
139+
xguide = "μ",
140+
yguide = "x")
141+
plot!(br_po, xguide = "μ", yguide = "x", label = "Maximum of periodic orbit")
142+
```
143+
144+
Let's see how to plot the periodic solution we just computed:
145+
146+
```@example Bif2
147+
sol = get_periodic_orbit(br_po, 10)
148+
plot(sol.t, sol[1,:], yguide = "x", xguide = "time", label = "")
149+
```
150+

0 commit comments

Comments
 (0)