From equilibria to equilibria

• From equilibria to equilibria
• • From simple branch point to equilibria
  • • Algorithm
    • Simple example: transcritical bifurcation
    • Simple example: pitchfork bifurcation
    • Simple example: quadratic parameter dependency
  • From non simple branch point to equilibria
  • • Examples
  • Assisted branching from non-simple bifurcation point
  • predictor(s)

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.

Algorithm

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

Simple example: transcritical bifurcation

using BifurcationKit, Plots

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

# parameters of the vector field
par = [-0.2]

# problem (with automatic differentiation)
prob = BifurcationProblem(F, [0.], par, 1; 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 = ["trivial branch", "bifurcated branch"], legend = :topleft)

Simple example: pitchfork bifurcation

using BifurcationKit, Plots

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

# parameters of the vector field
par = [-0.2]

# problem (with automatic differentiation)
prob = BifurcationProblem(F, [0.], par, 1; 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 = ["trivial branch", "bifurcated branch"], legend = :topleft)

Simple example: quadratic parameter dependency

In the previous examples, we need to have $O(p)$ dependency on the parameter $p$. If this is not the case, BifurcationKit can handle $O(p^2)$ dependency on the parameter $p$ although this is not generic. For higher order degeneracy, the user needs to provide the predictor him/herself.

using BifurcationKit, Plots

# vector field with degenerate transcritical bifurcation
Fdegenerate(x, p) = [x[1]^2 - p[1]^2]

# parameters of the vector field
par = [-1.0]

# problem (with automatic differentiation)
prob = BifurcationProblem(Fdegenerate, [1.], par, 1; record_from_solution = (x, p; k...)->x[1])

# compute branch of trivial equilibria
br = continuation(prob, PALC(), ContinuationPar(n_inversion = 6); normC = norminf)

# perform branch switching on both sides of the bifurcation point

br1 = continuation(br, 1, bothside = true)
scene = plot(br, br1; branchlabel = ["trivial branch", "bifurcated branch"], legend = :topleft)

Just for comparison, we can compute the reduced equation

bp = get_normal_form(br, 1)

Transcritical bifurcation point at p1 ≈ -2.8189256484623115e-17
Normal form (aδμ + b1⋅x⋅δμ + b2⋅x²/2 + b3⋅x³/6)
┌─ a  = 5.637851296924623e-17
├─ b1 = 0.0
├─ b2 = 2.0
└─ b3 = 0.0

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, LinearAlgebra

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.6, b = 1.5, c = 2.9)

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

# continuation options
opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.02, ds = 0.001,
	# 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 = 8)

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

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

- #  1,       nd at μ ≈ +0.00000281 ∈ (-0.00000065, +0.00000281), |δp|=3e-06, [converged], δ = ( 3,  0), step =  13
- #  2, endpoint at μ ≈ +0.40000000,                                                                     step =  28

You can now branch from the nd point

br2 = continuation(br, 1, opts_br; δp = 0.01)
plot(br, br2..., vars = (:param, :n))

For a better visualization, we can use the vector w defined above.

plot(br, br2..., vars = (:param, :nw))

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 μ ≈ 2.8052431305863613e-6. 
Kernel dimension = 3
Normal form:
 + 1.0 * x1 ⋅ p + -1.5 ⋅ x1³ + -2.9001 ⋅ x1 ⋅ x2² + -2.9001 ⋅ x1 ⋅ x3² + 0.6 ⋅ x2 ⋅ x3
 + 1.0 * x2 ⋅ p + -2.9001 ⋅ x1² ⋅ x2 + 0.6 ⋅ x1 ⋅ x3 + -1.5 ⋅ x2³ + -2.9001 ⋅ x2 ⋅ x3²
 + 1.0 * x3 ⋅ p + -2.9001 ⋅ x1² ⋅ x3 + 0.6 ⋅ 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.01
pred = predictor(bp, δp)

(before = [[0.0, 0.0, 0.0], [-0.05895671891836387, 0.05895671891836348, -0.058956718918364154]], after = [[0.0, 0.0, 0.0], [0.09640172394023658, 0.09640172394023658, 0.09640172394023658], [-0.09640172394023658, -0.09640172394023658, 0.09640172394023658], [-4.70197740328915e-37, -4.70197740328915e-37, 0.08164965809277261], [1.0965166582572064e-28, -1.1339875512552045e-28, -0.08164965809277261], [-1.4693679385278594e-39, -0.08164965809277261, 0.0]])

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 polynomial reduced equation bp(Val(:reducedForm), z, p).

These zeros serve as initial guesses for Newton's method applied to the full original functional:

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

(before = ┌─ Deflation operator with 2 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.009997194756869414, bpp = 0.010002805243130587)

These points can then be used to continue the different branches:

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

plot(br, brbp...)

predictor(s)

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