From equilibria to equilibria

• From equilibria to equilibria
• • From simple branch point to equilibria
  • • Simple example: transcritical bifurcation
    • Simple example: pitchfork bifurcation
    • Algorithm
  • From non simple branch point to equilibria
  • • Examples
  • Assisted branching from non-simple bifurcation point
  • predictors

From simple branch point to equilibria

See Branch switching (branch point) for the precise method definition

You can perform automatic branch switching by calling continuation with the following options:

continuation(br::ContResult, ind_bif::Int, optionsCont::ContinuationPar; kwargs...)

where br is a branch computed after a call to continuation with detection of bifurcation points enabled. This call computes the branch bifurcating from the ind_bifth bifurcation point in br. An example of use is provided in 2d Bratu–Gelfand problem.

See Branch switching (branch point) precise method definition

Simple example: transcritical bifurcation

using BifurcationKit, Plots

# vector field of transcritical bifurcation
F(x, p) = [x[1] * (p.μ - x[1])]

# parameters of the vector field
par = (μ = -0.2, )

# problem (automatic differentiation)
prob = BifurcationProblem(F, [0.], par, (@optic _.μ); record_from_solution = (x, p; k...) -> x[1])

# compute branch of trivial equilibria and detect a bifurcation point
br = continuation(prob, PALC(), ContinuationPar())

# perform branch switching on both sides of the bifurcation point
br1 = continuation(br, 1; bothside = true )

scene = plot(br, br1; branchlabel = ["br", "br1"], legend = :topleft)

Simple example: pitchfork bifurcation

using BifurcationKit, Plots

# vector field of pitchfork bifurcation
F(x, p) = [x[1] * (p.μ - x[1]^2)]

# parameters of the vector field
par = (μ = -0.2, )

# problem (automatic differentiation)
prob = BifurcationProblem(F, [0.], par, (@optic _.μ); record_from_solution = (x, p; k...) -> x[1])

# compute branch of trivial equilibria and
# detect a bifurcation point with improved precision
br = continuation(prob, PALC(), ContinuationPar(n_inversion = 6))

# perform branch switching on both sides of the bifurcation point
br1 = continuation(br, 1; bothside = true )

scene = plot(br, br1; branchlabel = ["br", "br1"], legend = :topleft)

Algorithm

• the normal form is computed and non-trivial zeros are used to produce guesses for points on the bifurcated branch.

From non simple branch point to equilibria

We provide an automatic branch switching method. The method is to first compute the reduced equation (see Non-simple branch point) and its zeros. These zeros are then seeded as initial guess for continuation. Hence, you can perform automatic branch switching by calling continuation with the following options:

continuation(br::ContResult, 
	ind_bif::Int,
	optionsCont::ContinuationPar;
	kwargs...)

Examples

An example of use is provided in 2d Bratu–Gelfand problem. A much simpler example is given now. It is a bit artificial because the vector field is its own normal form at the bifurcation point located at 0.

using BifurcationKit, Plots

function FbpD6(x, p)
    return [p.μ * x[1] + (p.a * x[2] * x[3] - p.b * x[1]^3 - p.c*(x[2]^2 + x[3]^2) * x[1]),
           p.μ * x[2] + (p.a * x[1] * x[3] - p.b * x[2]^3 - p.c*(x[3]^2 + x[1]^2) * x[2]),
           p.μ * x[3] + (p.a * x[1] * x[2] - p.b * x[3]^3 - p.c*(x[2]^2 + x[1]^2) * x[3])]
end

# model parameters
pard6 = (μ = -0.2, a = 0.3, b = 1.5, c = 2.9)

# problem
prob = BifurcationProblem(FbpD6, zeros(3), pard6, (@optic _.μ);
	record_from_solution = (x, p; k...) -> (n = norminf(x)))

# continuation options
opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.02, ds = 0.01,
	# parameter interval
	p_max = 0.4, p_min = -0.2,
	nev = 3,
	newton_options = NewtonPar(tol = 1e-10, max_iterations = 20),
	max_steps = 1000,
	n_inversion = 6)

br = continuation(prob, PALC(), opts_br)

 ┌─ Curve type: EquilibriumCont
 ├─ Number of points: 24
 ├─ Type of vectors: Vector{Float64}
 ├─ Parameter μ starts at -0.2, ends at 0.4
 ├─ Algo: PALC
 └─ Special points:

- #  1,       nd at μ ≈ +0.00008912 ∈ (-0.00002136, +0.00008912), |δp|=1e-04, [converged], δ = ( 3,  0), step =   8
- #  2, endpoint at μ ≈ +0.40000000,                                                                     step =  23

You can now branch from the nd point

br2 = continuation(br, 1, opts_br; δp = 0.02)

plot(br, br2...)

Assisted branching from non-simple bifurcation point

It may happen that the general procedure fails. We thus expose the procedure multicontinuation in order to let the user tune it to its need.

The first step is to compute the reduced equation, say of the first bifurcation point in br.

bp = get_normal_form(br, 1)

Non simple bifurcation point at μ ≈ 8.912198888088785e-5. 
Kernel dimension = 3
Normal form:
 + 1.0 * x1 ⋅ p + -1.5 ⋅ x1³ + -2.9001 ⋅ x1 ⋅ x2² + -2.9001 ⋅ x1 ⋅ x3² + 0.3 ⋅ x2 ⋅ x3
 + 1.0 * x2 ⋅ p + -2.9001 ⋅ x1² ⋅ x2 + 0.3 ⋅ x1 ⋅ x3 + -1.5 ⋅ x2³ + -2.9001 ⋅ x2 ⋅ x3²
 + 1.0 * x3 ⋅ p + -2.9001 ⋅ x1² ⋅ x3 + 0.3 ⋅ x1 ⋅ x2 + -2.9001 ⋅ x2² ⋅ x3 + -1.5 ⋅ x3³

Next, we want to find the zeros of the reduced equation. This is usually achieved by calling the predictor

δp = 0.005
pred = predictor(bp, δp)

(before = [[0.0, 0.0, 0.0]], after = [[0.0, 0.0, 0.0], [0.05382180219586076, 0.05382180219586075, 0.05382180219586075], [0.057735026918947635, 1.741994205345247e-25, 1.787117049123889e-25], [-0.057735026918947635, 8.609577764811676e-42, -1.1479437019748901e-41], [-2.8698592549372254e-42, -0.057735026918947635, -5.3809861030072976e-43], [5.605193857299268e-45, 0.057735026918947635, 7.006492321624085e-45], [-3.363116314379561e-44, 3.923635700109488e-44, -0.057735026918947635]])

which returns zeros of bp before and after the bifurcation point. You could also use your preferred procedure from Roots.jl (or other) to find the zeros of the polynomials bp(Val(:reducedForm), z, p).

We can use these zeros to form guesses to apply Newton for the full functional:

pts = BifurcationKit.get_first_points_on_branch(br, bp, pred, opts_br; δp)

(before = ┌─ Deflation operator with 1 root(s)
├─ eltype   = Float64
├─ power    = 2
├─ α        = 1.0
├─ dist     = dot
└─ autodiff = true
, after = ┌─ Deflation operator with 6 root(s)
├─ eltype   = Float64
├─ power    = 2
├─ α        = 1.0
├─ dist     = dot
└─ autodiff = true
, bpm = -0.004910878011119112, bpp = 0.005089121988880888)

We can then use this to continue the different branches

brbp = BifurcationKit.multicontinuation(br, bp, pts.before, pts.after, opts_br)

plot(br, brbp...)

Note that you can chose another predictor which uses all vertices of the cube as initial guesses

pred = predictor(bp, Val(:exhaustive), δp)
pts = BifurcationKit.get_first_points_on_branch(br, bp, pred, opts_br; δp)

(before = ┌─ Deflation operator with 1 root(s)
├─ eltype   = Float64
├─ power    = 2
├─ α        = 1.0
├─ dist     = dot
└─ autodiff = true
, after = ┌─ Deflation operator with 10 root(s)
├─ eltype   = Float64
├─ power    = 2
├─ α        = 1.0
├─ dist     = dot
└─ autodiff = true
, bpm = -0.004910878011119112, bpp = 0.005089121988880888)

brbp = BifurcationKit.multicontinuation(br, bp, pts.before, pts.after, opts_br)

plot(br, brbp...)

predictors

predictor(
    bp,
    δp;
    verbose,
    ampfactor,
    nbfailures,
    maxiter,
    perturb,
    J,
    normN,
    optn
)

predictor(
    bp,
    ,
    δp;
    verbose,
    ampfactor,
    nbfailures,
    maxiter,
    perturb,
    J,
    igs,
    normN,
    optn
)