# 🟡 Automatic diagram of 2d Bratu–Gelfand problem

The following example is exposed in Farrell, Patrick E., Casper H. L. Beentjes, and Ásgeir Birkisson. **The Computation of Disconnected Bifurcation Diagrams.** ArXiv:1603.00809 [Math], March 2, 2016. . It is also treated in Michiel Wouters. **Automatic Exploration Techniques for the Numerical Continuation of Large–Scale Nonlinear Systems**, 2019.

We consider the problem of Mittelmann:

\[\Delta u +NL(\lambda,u) = 0\]

with Neumann boundary condition on $\Omega = (0,1)^2$ and where $NL(\lambda,u)\equiv-10(u-\lambda e^u)$. This is a good example to show how automatic bifurcation diagram computation works.

We start with some imports:

```
using Revise, ForwardDiff
using BifurcationKit, LinearAlgebra, Plots, SparseArrays
const BK = BifurcationKit
# define the sup norm
norm2(x) = norm(x) / sqrt(length(x))
normbratu(x) = norm(x .* w) / sqrt(length(x)) # the weight w is defined below
# some plotting functions to simplify our life
plotsol!(x, nx = Nx, ny = Ny; kwargs...) = heatmap!(reshape(x, nx, ny); color = :viridis, kwargs...)
plotsol(x, nx = Nx, ny = Ny; kwargs...) = (plot();plotsol!(x, nx, ny; kwargs...))
```

and with the discretization of the problem

```
function Laplacian2D(Nx, Ny, lx, ly)
hx = 2lx/Nx
hy = 2ly/Ny
D2x = spdiagm(0 => -2ones(Nx), 1 => ones(Nx-1), -1 => ones(Nx-1) ) / hx^2
D2y = spdiagm(0 => -2ones(Ny), 1 => ones(Ny-1), -1 => ones(Ny-1) ) / hy^2
D2x[1,1] = -1/hx^2
D2x[end,end] = -1/hx^2
D2y[1,1] = -1/hy^2
D2y[end,end] = -1/hy^2
D2xsp = sparse(D2x)
D2ysp = sparse(D2y)
A = kron(sparse(I, Ny, Ny), D2xsp) + kron(D2ysp, sparse(I, Nx, Nx))
return A
end
ϕ(u, λ) = -10(u-λ*exp(u))
dϕ(u, λ) = -10(1-λ*exp(u))
function NL!(dest, u, p)
(;λ) = p
dest .= ϕ.(u, λ)
return dest
end
NL(u, p) = NL!(similar(u), u, p)
function Fmit!(f, u, p)
mul!(f, p.Δ, u)
f .= f .+ NL(u, p)
return f
end
Fmit(u, p) = Fmit!(similar(u), u, p)
```

It will also prove useful to have the jacobian of our functional and the other derivatives:

```
function JFmit(x,p)
J = p.Δ
dg = dϕ.(x, p.λ)
return J + spdiagm(0 => dg)
end
```

We need to pass the parameters associated to this problem:

```
Nx = 30
Ny = 30
lx = 0.5
ly = 0.5
# weight for normbratu
const w = (lx .+ LinRange(-lx,lx,Nx)) * (LinRange(-ly,ly,Ny))' |> vec
w .-= minimum(w)
Δ = Laplacian2D(Nx, Ny, lx, ly)
# parameters associated with the PDE
par_mit = (λ = .01, Δ = Δ)
# initial condition
sol0 = 0*ones(Nx, Ny) |> vec
# Bifurcation Problem
prob = BifurcationProblem(Fmit, sol0, par_mit, (@optic _.λ),; J = JFmit,
record_from_solution = (x, p; k...) -> (x = normbratu(x), n2 = norm(x), n∞ = norminf(x)),
plot_solution = (x, p; k...) -> plotsol!(x ; k...))
```

To compute the eigenvalues, we opt for the solver in `KrylovKit.jl`

```
# eigensolver
eigls = EigKrylovKit(dim = 70)
# options for Newton solver
opt_newton = NewtonPar(tol = 1e-8, verbose = true, eigsolver = eigls, max_iterations = 20)
# options for continuation, we want to locate very precisely the
# bifurcation points, so we tune the bisection accordingly
opts_br = ContinuationPar(dsmin = 0.0001, dsmax = 0.04, ds = 0.005, p_max = 3.5, p_min = 0.01, detect_bifurcation = 3, nev = 50, plot_every_step = 10, newton_options = (@set opt_newton.verbose = false), max_steps = 251, tol_stability = 1e-6, n_inversion = 6, dsmin_bisection = 1e-7, max_bisection_steps = 25, tol_bisection_eigenvalue = 1e-19)
```

Note that we put the option `detect_bifurcation = 3`

to detect bifurcations precisely with a **bisection** method. Indeed, we need to locate these branch points precisely to be able to call automatic branch switching.

In order to have an output like Auto07p, we provide the finaliser (see arguments of `continuation`

)

```
function finSol(z, tau, step, br; k...)
if length(br.specialpoint)>0
if br.specialpoint[end].step == step
BK._show(stdout, br.specialpoint[end], step)
end
end
return true
end
```

## Automatic bifurcation diagram

In order to avoid spurious branch switching, we use a callback (see `continuation`

) to reject specific continuation steps where the jump in parameters is too large or when the residual is too large:

```
function cb(state; kwargs...)
_x = get(kwargs, :z0, nothing)
fromNewton = get(kwargs, :fromNewton, false)
if ~fromNewton
return (norm(_x.u - state.x) < 20.5 && abs(_x.p - state.p) < 0.05)
end
true
end
```

Finally, before calling the automatic `bifurcationdiagram`

, we need to provide a function to adjust the continuation parameters as function of the branching level (Note that this function can be constant).

```
function optionsCont(x,p,l; opt0 = opts_br)
if l == 1
return opt0
elseif l==2
return setproperties(opt0 ;detect_bifurcation = 3,ds = 0.001, a = 0.75)
else
return setproperties(opt0 ;detect_bifurcation = 3,ds = 0.00051, dsmax = 0.01)
end
end
```

We are then ready to compute the bifurcation diagram. If we choose a level 5 of recursion like

```
diagram = @time bifurcationdiagram(prob, PALC(),
# important argument: this is the maximal
# recursion level
5,
optionsCont;
verbosity = 0, plot = true,
callback_newton = cb,
usedeflation = true,
finalise_solution = finSol,
normC = norminf)
```

this gives using `plot(diagram; plotfold = false, putspecialptlegend=false, markersize=2, title = "#branches = $(size(diagram))")`

:

We can zoom in on the left part to get

Actually, this plot is misleading because of the symmetries. If we chose a weighted norm which breaks those symmetries and use it to print the solution, we get

```
plot(diagram; plotfold = false, putspecialptlegend=false, markersize=2,
title = "#branches = $(size(diagram))", vars = (:param, :nw))
```

We can make more sense of these spaghetti by only plotting the first two levels of recursion

`plot(diagram; level = (1, 2), plotfold = false, putspecialptlegend=false, markersize=2, vars = (:param, :nw))`

## Interactive exploration

We can see that the non-simple 2d branch points (magenta points) have produced non trivial branches. For example, we can look at the second bifurcation point (the first is the fold) which is composed of 8 branches

`plot(getBranchesFromBP(diagram, 2); plotfold = false, legend = false, vars = (:param, :nw))`

## Interactive computation

Let's say you have been cautious and did not launch a deep bifurcation diagram computation by using a small recursion level 2:

```
diagram = bifurcationdiagram(prob, PALC(),
# here the recursion level is
2,
optionsCont;
verbosity = 0, plot = true,
callback_newton = cb,
usedeflation = true,
finalise_solution = finSol,
normC = norminf)
```

You would end up with this diagram

How can we complete this diagram without recomputing it from scratch? It is easy! For example, let us complete the magenta branches as follow

```
bifurcationdiagram!(prob,
# this improves the first branch on the violet curve. Note that
# for symmetry reasons, the first bifurcation point
# has 8 branches
get_branch(diagram, (1,)), 6, optionsCont;
verbosity = 0, plot = true,
callback_newton = cb,
finalise_solution = finSol,
usedeflation = true,
normC = norminf)
```

This gives the following diagram. Using this call, you can pinpoint the particular location where to refine the diagram.