Hutchinson with Diffusion

• Hutchinson with Diffusion
• • Problem discretization
  • Bifurcation analysis
  • Performance improvements
  • References

Consider the Hutchinson equation with diffusion

\[\begin{aligned} & \frac{\partial u(t, x)}{\partial t}=d \frac{\partial^2 u(t, x)}{\partial x^2}-a u(t-1, x)[1+u(t, x)], \quad t>0, x \in(0, \pi), \\ & \frac{\partial u(t, x)}{\partial x}=0, x=0, \pi \end{aligned}\]

where $a>0,d>0$.

Problem discretization

We start by discretizing the above PDE based on finite differences.

using Revise, DDEBifurcationKit, Parameters, LinearAlgebra, Plots, SparseArrays
using BifurcationKit
const BK = BifurcationKit

function Hutchinson(u, ud, p)
   @unpack a,d,Δ = p
   d .* (Δ*u) .- a .* ud[1] .* (1 .+ u)
end

delaysF(par) = [1.]

Bifurcation analysis

We can now instantiate the model

# discretisation
Nx = 200; Lx = pi/2;
X = -Lx .+ 2Lx/Nx*(0:Nx-1) |> collect

function DiffOp(N, lx; order = 2)
	hx = 2lx/N
	Δ = spdiagm(0 => -2ones(N), 1 => ones(N-1), -1 => ones(N-1) )
	Δ[1,1]=-1; Δ[end,end]=-1
	Δ = Δ / hx^2
	return Δ
end

Δ = DiffOp(Nx, Lx)

We are now ready to compute the bifurcation of the trivial (constant in space) solution:

# bifurcation problem
pars = (a = 0.5, d = 1, τ = 1.0, Δ = Δ, N = Nx)
x0 = zeros(Nx)

prob = ConstantDDEBifProblem(Hutchinson, delaysF, x0, pars, (@lens _.a))

optn = NewtonPar(eigsolver = DDE_DefaultEig())
opts = ContinuationPar(p_max = 10., p_min = 0., newton_options = optn, ds = 0.01, detect_bifurcation = 3, nev = 5, dsmax = 0.2, n_inversion = 4)
br = continuation(prob, PALC(), opts; verbosity = 0, plot = false, normC = norminf)
br

 ┌─ Curve type: EquilibriumCont
 ├─ Number of points: 42
 ├─ Type of vectors: Vector{Float64}
 ├─ Parameter a starts at 0.5, ends at 10.0
 ├─ Algo: PALC
 └─ Special points:

If `br` is the name of the branch,
ind_ev = index of the bifurcating eigenvalue e.g. `br.eig[idx].eigenvals[ind_ev]`

- #  1,     hopf at a ≈ +1.57668056 ∈ (+1.56784173, +1.57668056), |δp|=9e-03, [converged], δ = ( 2,  2), step =  10, eigenelements in eig[ 11], ind_ev =   2
- #  2,     hopf at a ≈ +2.26389996 ∈ (+2.26169025, +2.26389996), |δp|=2e-03, [converged], δ = ( 2,  2), step =  13, eigenelements in eig[ 14], ind_ev =   4
- #  3,     hopf at a ≈ +4.75534651 ∈ (+4.75424166, +4.75534651), |δp|=1e-03, [converged], δ = ( 2,  2), step =  22, eigenelements in eig[ 23], ind_ev =   6
- #  4,     hopf at a ≈ +9.43550952 ∈ (+9.43109010, +9.43550952), |δp|=4e-03, [converged], δ = ( 2,  2), step =  39, eigenelements in eig[ 40], ind_ev =   8
- #  5, endpoint at a ≈ +10.00000000,                                                                      step =  41

We note that the first Hopf point is close to the theoretical value $a=\frac\pi 2$. This can be improved by increasing opts.n_inversion.

We can now plot the branch

scene = plot(br)

Performance improvements

The previous implementation being simple, it leaves a lot performance on the table. For example, the jacobian is dense because it is computed with automatic differentiation without sparsity detection.

We show how to specify the jacobian and speed up the code a lot.

# analytical jacobian
function JacHutchinson(u, p)
   @unpack a,d,Δ = p
   # we compute the jacobian at the steady state
   J0 = d * Δ .- a .* Diagonal(u)
   J1 = -a .* Diagonal(1 .+ u)
   return J0, [J1]
end

prob2 = ConstantDDEBifProblem(Hutchinson, delaysF, x0, pars, (@lens _.a); J = JacHutchinson)

optn = NewtonPar(eigsolver = DDE_DefaultEig())
opts = ContinuationPar(p_max = 10., p_min = 0., newton_options = optn, ds = 0.01, detect_bifurcation = 3, nev = 5, dsmax = 0.2, n_inversion = 4)
br = continuation(prob2, PALC(), opts; verbosity = 1, plot = true, normC = norminf)
br

 ┌─ Curve type: EquilibriumCont
 ├─ Number of points: 42
 ├─ Type of vectors: Vector{Float64}
 ├─ Parameter a starts at 0.5, ends at 10.0
 ├─ Algo: PALC
 └─ Special points:

If `br` is the name of the branch,
ind_ev = index of the bifurcating eigenvalue e.g. `br.eig[idx].eigenvals[ind_ev]`

- #  1,     hopf at a ≈ +1.57668056 ∈ (+1.56784173, +1.57668056), |δp|=9e-03, [converged], δ = ( 2,  2), step =  10, eigenelements in eig[ 11], ind_ev =   2
- #  2,     hopf at a ≈ +2.26389996 ∈ (+2.26169025, +2.26389996), |δp|=2e-03, [converged], δ = ( 2,  2), step =  13, eigenelements in eig[ 14], ind_ev =   4
- #  3,     hopf at a ≈ +4.75534651 ∈ (+4.75424166, +4.75534651), |δp|=1e-03, [converged], δ = ( 2,  2), step =  22, eigenelements in eig[ 23], ind_ev =   6
- #  4,     hopf at a ≈ +9.43550952 ∈ (+9.43109010, +9.43550952), |δp|=4e-03, [converged], δ = ( 2,  2), step =  39, eigenelements in eig[ 40], ind_ev =   8
- #  5, endpoint at a ≈ +10.00000000,                                                                      step =  41

References