BifurcationKit.jl: Automatic Bifurcation Diagrams in Julia

[rveltz]

Dear All,

I am writing this post to report about my package BifurcationKit.jl. First, this is a new name for an old package

The focus of this new version is automatic bifurcation diagram (aBD) computation for PDEs. Although this is a bit preliminary, I am excited to share some advances regarding the case of equilibria (see figure below which shows the equilibria of a PDE as function of a parameter (x-axis)).

Notably, the algorithm (aBD) is quite short, what it requires is a robust branch switching algorithm. As you can tell, computing this diagram manually (by telling which branch point to consider) will be quite time consuming and that’s why you usually find partial version of this in the litterature.

There aren’t so many softwares which allow to do this, and even less for PDE or large scale problems (I know one which is not open source nor released). For now, the aBD works for simple branch points but I will soon commit a version which works even when the PDE/functional has symmetries (see teaser). With this, the package would definitely be beyond state of the art…

You can find a basic example and a more involved one for a 1d PDE.

Hopefully, some of you will find this useful,

Bests

[daybrown]

Hi, This is a great package for bifurcation, thank you.
I am trying to make it work for a simple ODE system, but my hopf continuation failed. Below is my code:

function Fitz(u,p)
    x,y = u
    @unpack I,a,b,c = p
    du = similar(u)
    
    du[1] = x-x^3/3-y + I
    du[2] = a*(b*x-c*y)
    
    return du
end

sol0 = [0.1,0.1]

par_Fitz = (I=-2.0,a=0.08,b=1.,c=0.8)

ls = GMRESKrylovKit(dim = 100)
optnewton = NewtonPar(tol = 1e-11, verbose = true, linsolver = ls)
out, _, _ = @time newton( Fitz, sol0, par_Fitz, optnewton)

optnew = NewtonPar(verbose = true, tol = 1e-10)
optcont = ContinuationPar(dsmin = 0.01, dsmax = 0.1, ds= 0.05, pMax = 5.1, pMin = -2.1,
    newtonOptions = setproperties(optnew; maxIter = 50, tol = 1e-10), 
	maxSteps = 300, plotEveryStep = 40, 
	detectBifurcation = 3, nInversion = 4, tolBisectionEigenvalue = 1e-12, dsminBisection = 1e-5)

br, _ = @time continuation(Fitz, out, par_Fitz, (@lens _.I),
		optcont; plot = true, verbosity = 0,
		printSolution = (x, p) -> x[1])


Ddt(f, x, p, dx) = ForwardDiff.derivative(t -> f(x .+ t .* dx, p), 0.)
dFbpSecBif(x,p)         =  ForwardDiff.jacobian( z-> Fitz(z,p), x)
d1FbpSecBif(x,p,dx1)         = Ddt((z, p0) -> Fitz(z, p0), x, p, dx1)
d2FbpSecBif(x,p,dx1,dx2)     = Ddt((z, p0) -> d1FbpSecBif(z, p0, dx1), x, p, dx2)
d3FbpSecBif(x,p,dx1,dx2,dx3) = Ddt((z, p0) -> d2FbpSecBif(z, p0, dx1, dx2), x, p, dx3)
jet = (Fitz, dFbpSecBif, d2FbpSecBif, d3FbpSecBif)

# index of the Hopf point in br.bifpoint
norminf = x -> norm(x, Inf)
ind_hopf = 1
hopfpoint, _, flag = @time newton(Fitz, dFbpSecBif,
	br, ind_hopf, par_Fitz, (@lens _.I); normN = norminf)
flag && printstyled(color=:red, "--> We found a Hopf Point at μ = ", hopfpoint.p[1],
    ", from μ = ", br.bifpoint[ind_hopf].param, "\n")

# automatic branch switching from Hopf point
opt_po = NewtonPar(tol = 1e-5, verbose = true, maxIter = 15)
opts_po_cont = ContinuationPar(dsmin = 0.01, dsmax = 0.1, ds = 0.1, pMax = 5.5,pMin=-2.1, maxSteps = 200,
    newtonOptions = opt_po, saveSolEveryStep = 2, plotEveryStep = 1, nev = 11, precisionStability = 1e-4,
	detectBifurcation = 3, dsminBisection = 1e-4, maxBisectionSteps = 10);

# number of time slices for the periodic orbit
M = 51
probFD = PeriodicOrbitTrapProblem(M = M)
br_po, = continuation(
	# arguments for branch switching from the first 
	# Hopf bifurcation point
	jet..., br, 1,
	# arguments for continuation
	opts_po_cont, probFD;
	# OPTIONAL parameters
	# we want to jump on the new branch at phopf + δp
	# ampfactor is a factor to increase the amplitude of the guess
	δp = 0.05, ampfactor = 1.1,
	# specific method for solving linear system
	# of Periodic orbits with trapeze method
	linearPO = :FullLU,
# regular options for continuation
	verbosity = 3,	plot = true, 
	printSolution = (x,p)->x[1], normC = norminf)

The branch br is correct, while br_po is not correct. Could you tell what is wrong?

[rveltz]

Thank you for trying the package!

You should open an issue on BifurcationKit.jl, I will be able to help you more there.

First, you are using an iterative linear solver, use ls = defaultLS() instead.

You are far from the hopf point with your options, you should try:

usedeflation = true, # find nontrivial periodic solution
# OPTIONAL parameters
# we want to jump on the new branch at phopf + δp
# ampfactor is a factor to increase the amplitude of the guess
δp = 0.001, ampfactor = 1.0,

Also, the Hopf normal form shows that the period is about 22, you should try M = 150. The problem seems quite stiff, and the period will grow a lot, making the computation of periodic orbits with trapezoidal method quite difficult ; we dont have adaptive time step yet. An alternative is to try shooting methods then. I get this with shooting:

Additionally, you are looking at a notoriously difficult problem here with canard explosion. The methods in BifurcationKit.jl for periodic orbits, (for ODE), are not in par with established software like Auto or MatCont which use adaptive collocation methods. However these methods do not scale well for PDE, so I did not implement them. You can look at GitHub - tkf/Bifurcations.jl: Bifurcation analysis in pure Julia. Totally work-in-progress. also

[daybrown]

Hi,

Thank you for the reply, and I saw the result you got, which looks good.
I opened an issue in BifurcationKit.jl, could you help with modifying my code to get a similar result?