# Automatic bifurcation diagram for a simple model

**URL:** https://discourse.julialang.org/t/automatic-bifurcation-diagram-for-a-simple-model/87227
**Category:** Numerics
**Created:** [September 13, 2022, 5:43pm UTC](https://discourse.julialang.org/t/automatic-bifurcation-diagram-for-a-simple-model/87227 "2022-09-13T17:43:11Z")
**Posts on this page:** 1
**Page:** 1

<div class="post-metadata">

### Author: ![zmoitier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/zmoitier/32/29102_2.png) [@zmoitier](https://discourse.julialang.org/u/zmoitier)
#### Post date: [September 13, 2022, 5:43pm UTC](https://discourse.julialang.org/t/automatic-bifurcation-diagram-for-a-simple-model/87227/1 "2022-09-13T17:43:11Z")

</div>

I am trying to compute the bifurcation diagram using `BifucationKit.jl` associated to the equation: find a vector `u` and a scalar `λ` such that `Au - λu - u.^3 = 0` where `A` is a real symmetric matrix. For this problem, the theory tell us that on the trivial branch `(u=0, λ)` only the eigenvalues of the matrix `A` are bifurcation points. When I adapt the code found in [example](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/BifurcationDiagram/) to this problem I get spurious bifurcation points.

```julia
using LinearAlgebra, Setfield, Parameters
using BifurcationKit
using Plots
const BK = BifurcationKit

# discretisation
N = 32
A = ((N + 1) / π)^2 * SymTridiagonal(fill(2, N), fill(-1, N - 1))

display(eigvals(A)[1:3])

# The function F and its differential
function F(u::Vector, param::NamedTuple)
    return param.A * u - param.λ .* u - u .^ 3
end

function DF(u::Vector, param::NamedTuple)
    return param.A - UniformScaling(param.λ) - Diagonal(3 .* u .^ 2)
end

# parameters associated with the equation
param = (λ=0.0, A=A)
λmin, λmax = -5.0, 5.0

# bifurcation problem
prob = BifurcationProblem(F, zeros(N), param,
                          # specify the continuation parameter
                          (@lens _.λ); J=DF, recordFromSolution=(u, p) -> norm(u))

# options for Krylov-Newton
opt_newton = NewtonPar(; tol=1e-9, maxIter=20)

# options for continuation
opts_br = ContinuationPar(; dsmin=0.001, dsmax=0.05, ds=0.01, θ=0.5,
                          # parameter interval
                          pMin=λmin, pMax=λmax,
                          maxSteps=200, newtonOptions=opt_newton,
                          nev=4,
                          # detect bifurcations with bisection method
                          detectBifurcation=3, nInversion=4, tolBisectionEigenvalue=1e-8,
                          dsminBisection=1e-9)

diagram = bifurcationdiagram(prob, PALC(), 2,
                             (args...) -> setproperties(opts_br; pMin=λmin, pMax=λmax,
                                                        ds=0.001, dsmax=0.005, nInversion=8,
                                                        detectBifurcation=3,
                                                        dsminBisection=1e-18,
                                                        maxBisectionSteps=20))

# You can plot the diagram like
plot(diagram; putspecialptlegend=false, markersize=2, plotfold=false,
     title="#branches = $(size(diagram))")

display(diagram)

```

It found 3 bifurcation points, the first two should not exists and the third is a true bifurcation point but the new branch is not computed on the diagram. I must be doing something wrong but I do know what.

```julia
────────────────────────────────────────────────────────────────────────────────
--> New branch, level = 2, dim(Kernel) = 3, code = 0, from bp #1 at p = 0.001627588557867863, type = nd
- # 1, nd at p ≈ +0.00162759 ∈ (+0.00162759, +0.00162759), |δp|=3e-09, [guessL], δ = (-3, 0), step = 2, eigenelements in eig[3], ind_ev = 32
────────────────────────────────────────────────────────────────────────────────
--> New branch, level = 2, dim(Kernel) = 3, code = 0, from bp #2 at p = 3.6602852057769955e-5, type = nd
- # 2, nd at p ≈ +0.00003660 ∈ (+0.00003660, +0.00162759), |δp|=2e-03, [guess], δ = ( 3, 0), step = 3, eigenelements in eig[4], ind_ev = 32
────────────────────────────────────────────────────────────────────────────────
--> New branch, level = 2, dim(Kernel) = 1, code = 0, from bp #3 at p = 0.9992458117303361, type = bp
- # 3, bp at p ≈ +0.99924581 ∈ (+0.99924409, +0.99924581), |δp|=2e-06, [converged], δ = (-1, 0), step = 145, eigenelements in eig[146], ind_ev = 32
┌ Warning: The bifurcating eigenvalue is not that close to Re = 0. We found 424.54824195236205 !≈ 0. You can perhaps increase the argument `nev`.
└ @ BifurcationKit ~/.julia/packages/BifurcationKit/fnCOr/src/NormalForms.jl:33

```

* * *

Additional information:

- Julia: version `v1.8.1`
- `BifurcationKit.jl`: version `v0.2.0`