# Getting Started with BifurcationKit

This tutorial will introduce you to the functionalities for computing bifurcation diagrams and follow branches of solutions.

## Example 1: solving the perturbed pitchfork equation

In this example, we will solve the equation

\[\mu + x-\frac{x^3}{3}=0\]

as function of $\mu$ by looking at the solutions in the connected component of $(x_0,\mu_0)\approx(-2,-1)$. Here $x\in\mathbb R$ is the state variable and $\mu$ is our parameter. The general workflow is to define a problem, solve the problem, and then analyze the solution. The full code for solving this problem is:

```
using BifurcationKit, Plots
F(x, p) = @. p.μ + x - x^3/3
prob = BifurcationProblem(F, [-2.], (μ = -1.,), (@lens _.μ);
record_from_solution = (x,p) -> (x = x[1]))
br = continuation(prob, PALC(), ContinuationPar(p_min = -1., p_max = 1.))
plot(br)
```

where the pieces are described below.

### Step 1: Defining a problem

To solve this numerically, we define a problem type by giving it the equation, the initial condition, the parameters and the parameter axis to solve over:

```
using BifurcationKit
F(x, p) = @. p.μ + x - x^3/3
prob = BifurcationProblem(F, [-2.], (μ = -1.,), (@lens _.μ);
record_from_solution = (x,p) -> (x = x[1]))
```

```
┌─ Bifurcation Problem with uType Vector{Float64}
├─ Inplace: false
├─ Symmetric: false
└─ Parameter: μ
```

Note that BifurcationKit.jl will choose the types for the problem based on the types used to define the problem type. For our example, notice that `u0`

is a `Vector{Float64}`

, and therefore this will solve with the dependent variables being `Vector{Float64}`

. You can use this to choose to solve with `Float32`

for example to run this on the GPU (see example).

You can customize a few scalar indicators for each step (for example if you don't want to save all solutions) by providing a function `record_from_solution`

. You can also control how the solution is plotted during a continuation run by providing a function `plot_solution`

. This is especially useful when studying PDE for example.

### Step 2: Solving a problem

After defining a problem, you "solve" it using `continuation`

.

`br = continuation(prob, PALC(), ContinuationPar(p_min = -1., p_max = 1.))`

The solvers can be controlled using the available options `ContinuationPar`

. For example, we can increase the maximum continuation step (in order to get a less points) by using the command `dsmax = 0.25`

```
using Plots
opts = ContinuationPar(p_min = -1., p_max = 1., dsmax = 0.25, max_steps = 1000)
br = continuation(prob, PALC(), opts)
scene = plot(br)
```

### Choosing a continuation algorithm

BifurcationKit.jl offers a much wider variety of continuation algorithms than traditional continuation softwares. Many of these algorithms are from recent research and have their own strengths and weaknesses. Each algoritm comes with a doc string, for example:

`BifurcationKit.PALC`

— Type`struct PALC{Ttang<:BifurcationKit.AbstractTangentComputation, Tbls<:BifurcationKit.AbstractLinearSolver, T, Tdot} <: BifurcationKit.AbstractContinuationAlgorithm`

Pseudo-arclength continuation algorithm.

Additional information is available on the website.

**Fields**

`tangent::BifurcationKit.AbstractTangentComputation`

: Tangent predictor, must be a subtype of`AbstractTangentComputation`

. For example`Secant()`

or`Bordered()`

, Default: Secant()`θ::Any`

:`θ`

is a parameter in the arclength constraint. It is very**important**to tune it. It should be tuned for the continuation to work properly especially in the case of large problems where the < x - x*0, dx*0 > component in the constraint equation might be favoured too much. Also, large thetas favour p as the corresponding term in N involves the term 1-theta. Default: 0.5`_bothside::Bool`

: [internal], Default: false`bls::BifurcationKit.AbstractLinearSolver`

: Bordered linear solver used to invert the jacobian of the newton bordered problem. It is also used to compute the tangent for the predictor`Bordered()`

, Default: MatrixBLS()`dotθ::Any`

:`dotθ = DotTheta()`

, this sets up a dot product`(x, y) -> dot(x, y) / length(x)`

used to define the weighted dot product (resp. norm) $\|(x, p)\|^2_\theta$ in the constraint $N(x, p)$ (see online docs on PALC). This argument can be used to remove the factor`1/length(x)`

for example in problems where the dimension of the state space changes (mesh adaptation, ...) Default: DotTheta()

For example, you can chose a different tangent predictor in `PALC`

```
opts = ContinuationPar(p_min = -1., p_max = 1.)
br = continuation(prob, PALC(tangent = Bordered()), opts)
scene = plot(br)
```

or you can use the Moore-Penrose continuation algorithm

```
opts = ContinuationPar(p_min = -1., p_max = 1.)
br = continuation(prob, MoorePenrose(), opts)
scene = plot(br)
```

### Step 3: Analyzing the solution

The result of `continuation`

is a solution object. A summary of the result is provided by the `show`

method:

`show(br) # this is equivalent to the REPL julia> br`

```
┌─ Curve type: EquilibriumCont
├─ Number of points: 55
├─ Type of vectors: Vector{Float64}
├─ Parameter μ starts at -1.0, ends at 1.0
├─ Algo: MoorePenrose
└─ Special points:
- # 1, bp at μ ≈ +0.66597167 ∈ (+0.66597167, +0.66658738), |δp|=6e-04, [converged], δ = ( 1, 0), step = 20
- # 2, bp at μ ≈ -0.66665277 ∈ (-0.66666619, -0.66665277), |δp|=1e-05, [converged], δ = (-1, 0), step = 38
- # 3, endpoint at μ ≈ +1.00000000, step = 54
```

From there, you can see that the branch has 55 points, the algorithm is also recalled because it can be modified internally. This summary shows that two bifurcation points where detected. At each such point, the couple `δ`

indicates how many real/complex eigenvalues crossed the imaginary axis. This is useful for debugging or when non generic bifurcations are encountered.

We can access the 5th value of the branch with:

`br[5]`

`(x = -2.0742615684285792, param = -0.9006174451276318, itnewton = 2, itlinear = 0, ds = 0.04102643120243098, n_unstable = 0, n_imag = 0, stable = true, step = 4, eigenvals = ComplexF64[-3.3025610542597894 + 0.0im], eigenvecs = ComplexF64[1.0 + 0.0im;;])`

The solution contains many other fields:

`propertynames(br)`

`(:branch, :eig, :sol, :contparams, :kind, :prob, :specialpoint, :alg)`

Hence, the eigenelements are saved in `br.eig`

, the solutions are saved in `br.sol`

and the bifurcation points in `br.specialpoint`

.

### Plotting branches

While one can directly plot solution time points using the tools given above, convenience commands are defined by recipes for Plots.jl. To plot the solution object, simply call plot:

```
#]add Plots # You need to install Plots.jl before your first time using it!
using Plots
#plotly() # You can optionally choose a plotting backend
plot(br)
```

## Example 2: simple branching

In this example, we will solve the equation

\[0 = x\cdot(\mu-x)\]

as function of $\mu$. Here $u\in\mathbb R$ is the state variable and $\mu$ is our parameter.

In our example, we know by calculus that the solutions to this equation are $u_0(\mu)=0$ and $u_1(\mu)=\mu$ but we will use BifurcationKit.jl to solve this problem numerically, which is essential for problems where a symbolic solution is not known.

In case we know there are many branches, the best is to use an automatic method to compute them all. We will focus on `bifurcationdiagram`

which computes the connected component of the initial guess in the plane $(x,\mu)$. An alternative is to use Deflated Continuation.

We define a problem type by giving it the equation, the initial condition, the parameters and the parameter axis to solve over:

```
using Plots
using BifurcationKit
Fbp(u, p) = @. u * (p.μ - u)
# bifurcation problem
prob = BifurcationProblem(Fbp, [0.0], (μ = -0.2,),
# specify the continuation parameter
(@lens _.μ),
record_from_solution = (x, p) -> x[1])
```

```
┌─ Bifurcation Problem with uType Vector{Float64}
├─ Inplace: false
├─ Symmetric: false
└─ Parameter: μ
```

We then aim at calling `bifurcationdiagram`

which will do the jobs of computing recursively the branches which are connected together. Compared to `continuation`

, `bifurcationdiagram`

requires the maximal level of recursion (in this case 2 because there are 2 branches) and a function providing the continuation parameters for each branch (which may differ from branch to branch if the user decides). This explains the following code:

```
# options for continuation
opts_br = ContinuationPar(
# parameter interval
p_max = 0.2, p_min = -0.2,
# detect bifurcations with bisection method
# we increase the precision of the bisection
n_inversion = 4)
# automatic bifurcation diagram computation
diagram = bifurcationdiagram(prob, PALC(),
# very important parameter. This specifies the maximum amount of recursion
# when computing the bifurcation diagram. It means we allow computing branches of branches
# at most in the present case.
2,
(args...) -> opts_br,
)
```

```
[Bifurcation diagram]
┌─ From 0-th bifurcation point.
├─ Children number: 2
└─ Root (recursion level 1)
┌─ Curve type: EquilibriumCont
├─ Number of points: 9
├─ Type of vectors: Vector{Float64}
├─ Parameter μ starts at -0.2, ends at 0.2
├─ Algo: PALC
└─ Special points:
- # 1, bp at μ ≈ +0.00108349 ∈ (-0.00333593, +0.00108349), |δp|=4e-03, [converged], δ = ( 1, 0), step = 6
- # 2, endpoint at μ ≈ +0.20000000, step = 8
```

You can plot the diagram like

`plot(diagram)`

## Example 3: continuing periodic orbits

In this example, we will compute periodic orbits of the Stuart-Landau oscillator:

\[\begin{aligned} \frac{du}{dt} &= r u - \nu v - (u^2+v^2) (c_3 u - \mu v) \\ \frac{dv}{dt} &= r v + \nu u - (u^2+v^2) (c_3 + \mu u). \end{aligned}\]

The ODE is easily written with a function:

```
using BifurcationKit, Plots
function Fsl(X, p)
(;r, μ, ν, c3) = p
u, v = X
ua = u^2 + v^2
[
r * u - ν * v - ua * (c3 * u - μ * v)
r * v + ν * u - ua * (c3 * v + μ * u)
]
end
```

`Fsl (generic function with 1 method)`

and then we can use this to define a bifurcation problem:

```
par_sl = (r = 0.1, μ = 0., ν = 1.0, c3 = 1.0)
u0 = zeros(2)
prob = BifurcationProblem(Fsl, u0, par_sl, (@lens _.r))
```

```
┌─ Bifurcation Problem with uType Vector{Float64}
├─ Inplace: false
├─ Symmetric: false
└─ Parameter: r
```

For this simple problem, we detect the existence of periodic orbits by locating a Hopf bifurcation. This is done as in the previous example by continuing the zero solution:

`br = continuation(prob, PALC(), ContinuationPar(), bothside = true)`

```
┌─ Curve type: EquilibriumCont
├─ Number of points: 25
├─ Type of vectors: Vector{Float64}
├─ Parameter r starts at -1.0, ends at 1.0
├─ Algo: PALC
└─ Special points:
- # 1, endpoint at r ≈ -1.00000000, step = 0
- # 2, hopf at r ≈ -0.00595553 ∈ (-0.00595553, +0.00299379), |δp|=9e-03, [converged], δ = ( 2, 2), step = 8
- # 3, endpoint at r ≈ +1.00000000, step = 24
```

In the result above, we see that a Hopf bifurcation has been detected:

`- # 2, hopf at r ≈ -0.00595553 ∈ (-0.00595553,...`

We compute the branch of periodic orbits which is nearby. We thus provide the branch `br`

, the index of the special point we want to branch from: 2 in this case and a method `PeriodicOrbitOCollProblem(20, 5)`

to compute periodic orbits. You can look at Periodic orbits computation for a list of all methods. Suffice it to say that `PeriodicOrbitOCollProblem`

is the default method in the case of ODEs.

```
br_po = continuation(br, 2, ContinuationPar(),
PeriodicOrbitOCollProblem(20, 5)
)
```

```
┌─ Curve type: PeriodicOrbitCont from Hopf bifurcation point.
├─ Number of points: 15
├─ Type of vectors: Vector{Float64}
├─ Parameter r starts at 0.00404446825657828, ends at 1.0
├─ Algo: PALC
└─ Special points:
- # 1, endpoint at r ≈ +1.00000000, step = 14
```

### Analyzing the solution

The branch of periodic orbits has been computed. You can look at what is recorded in the first point on the branch:

`br_po[1]`

`(max = 0.06359613398726449, min = -0.0635961339872645, amplitude = 0.12719226797452898, period = 6.283185307179591, param = 0.00404446825657828, itnewton = 0, itlinear = 0, ds = 0.01, n_unstable = 0, n_imag = 0, stable = true, step = 0, eigenvals = ComplexF64[-9.73221503391148e-12 + 0.0im, -0.05082428707935699 + 0.0im], eigenvecs = ComplexF64[-3.920726923383661e-14 + 0.0im -1.0 + 0.0im; -1.0 + 0.0im -6.721246154371993e-15 + 0.0im])`

It shows that the maximum/minimum/amplitude/period of the periodic orbit are recorded by default. You can also plot all the branches as follows

`plot(br, br_po, branchlabel = ["equilibria", "periodic orbits"])`

Finally, if you are interested in the periodic orbits saved in `br_po`

, for example to plot it, the method `get_periodic_orbit`

is what you are looking for:

```
sol = get_periodic_orbit(br_po, 10)
plot(sol.t, sol[1,:], label = "u", xlabel = "time")
plot!(sol.t, sol[2,:], label = "v", xlabel = "time")
```

### Plotting the periodic orbit during continuation

If you plot the solution during continuation, you see that the right bottom panel is empty ; this panel is used to plot the solution at the current continuation step:

```
br_po = continuation(br, 2, opts,
PeriodicOrbitOCollProblem(20, 5);
plot = true,
)
scene = title!("") #hide
```

(Note that the bottom panel is a plot of the eigenvalues of the jacobian in the complex plane at the current continuation step. )

In order to plot the periodic solution during continuation, you need to supply a `periodic_solution`

to `continuation`

. This is not done by default because in some cases, obtaining the solution is costly (*e.g.* for Shooting methods). Based on the previous paragraph, it is straightforward to implement this plotting function:

```
br_po = continuation(br, 2, opts,
PeriodicOrbitOCollProblem(20, 5);
plot = true,
plot_solution = (x, par; k...) -> begin
# par is a Named tuple which contains
# the problem for computing periodic orbits
# and the value of the parameter at the current step
sol = get_periodic_orbit(par.prob, x, par.p)
plot!(sol.t, sol.u'; xlabel = "time", label="", k...)
end
)
scene = title!("") #hide
```