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 to this problem I get spurious bifurcation points.
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.
────────────────────────────────────────────────────────────────────────────────
--> 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: versionv0.2.0