🟠 Fronts in 1d autocatalytic model (Automatic)

• 🟠 Fronts in 1d autocatalytic model (Automatic)
• • Freezing method
  • Continuation of front solutions
  • References

We consider the following model ^{[Balmforth][Malham]} which is also treated in ^[Beyn]

\[\begin{array}{l} u_{t}=a u_{x x}-u f(v), \quad a>0, u, v: \mathbb{R} \rightarrow \mathbb{R} \\ v_{t}=v_{x x}+u f(v) \end{array}\]

where $f(u) = u^m 1_{u\geq 0}$. We chose the boundary conditions

\[\left(u_{-}, v_{-}\right)=(0,1),\quad \left(u_{+}, v_{+}\right)=(1,0)\tag{BC}.\]

It is straightforward to implement this problem as follows:

using Revise
using ForwardDiff, SparseArrays
using BifurcationKit, LinearAlgebra, Plots
const BK = BifurcationKit

# supremum norm
f(u) = u^9 # solutions are positive, so remove the heaviside

# helper function to plot solutions
function plotsol!(x; k...)
	u = @view x[1:end÷2]
	v = @view x[end÷2:end]
	plot!(u; label="u", k...)
	plot!(v; label="v", k...)
end
plotsol(x; k...) = (plot();plotsol!(x; k...))

# encode the nonlinearity
@views function NL!(dest, U, p, t = 0.)
	N = p.N
	u = U[1:N]
	v = U[N+1:2N]
	dest[1:N]    .= -u .* f.(v)
	dest[N+1:2N] .=  -dest[1:N]#u .* f.(v)
	return dest
end

# function for the differential with specific boundary conditions
# for fronts
@views function applyD_add!(f, U, p, a)
	uL = 0; uR = 1; vL = 1; vR = 0
	n = p.N
	u = U[1:n]
	v = U[n+1:2n]

	c1 = 1 / (2p.h)
	f[1]   += a * (uL      - u[2] ) * c1
	f[end] += a * (v[n-1]  - vR   ) * c1

	f[n]   += a * (u[n-1] - uR  ) * c1
	f[n+1] += a * (    vL - v[2] ) * c1

	@inbounds for i=2:n-1
		  f[i] += a * (u[i-1] - u[i+1]) * c1
		f[n+i] += a * (v[i-1] - v[i+1]) * c1
	end
	return f
end

# function which encodes the full PDE
@views function Fcat!(f, U, p, t = 0)
	uL = 0; uR = 1; vL = 1; vR = 0
	n = p.N
	# nonlinearity
	NL!(f, U, p)

	# Dirichlet boundary conditions
	h2 = p.h * p.h
	c1 = 1 / h2

	u = U[1:n]
	v = U[n+1:2n]

	f[1]   += p.a * (uL      - 2u[1] + u[2] ) * c1
	f[end] +=       (v[n-1]  - 2v[n] + vR   ) * c1

	f[n]   += p.a * (u[n-1] - 2u[n] +  uR  ) * c1
	f[n+1] +=       (vL - 2v[1] + v[2] ) * c1

	@inbounds for i=2:n-1
		  f[i] += p.a * (u[i-1] - 2u[i] + u[i+1]) * c1
		f[n+i] +=       (v[i-1] - 2v[i] + v[i+1]) * c1
	end
	return f
end
Fcat(x, p, t = 0.) = Fcat!(similar(x), x, p, t)
Jcat(x,p) = sparse(ForwardDiff.jacobian(x -> Fcat(x, p), x))

We chose the following parameters:

N = 200
lx = 25.
X = LinRange(-lx,lx, N)
par_cat = (N = N, a = 0.18, h = 2lx/N)

u0 = @. (tanh(2X)+1)/2
U0 = vcat(u0, 1 .- u0)

# we define a problem to hold the vector field
prob = BifurcationProblem(Fcat, u0, par_cat, (@optic _.a); J = Jcat)

Freezing method

The problem may feature fronts, that is solutions of the form $u(x,t) = \tilde u(x-st)$ (same for $v$) for a fixed value of the profile $\tilde u$ and the speed $s$. The equation for the front profile is, up to an abuse of notations (we removed the tildes)

\[\begin{array}{l} 0=a u_{\xi\xi}+s\cdot u_{\xi}-u f(v)\\ 0=v_{\xi\xi}+s\cdot v_{\xi}+u f(v) \end{array}\]

with unknowns $u,v,s$. The front is solution of these equations but it is not uniquely determined because of the phase invariance. Hence, we add the phase condition (see ^[Beyn])

\[0 = \left\langle (u,v), \partial_\xi (u_0,v_0) \right\rangle\]

where $U_0:=(u_0,v_0)$ is some fixed profile. This is encoded in the problem TWProblem

using LinearAlgebra

# this structure encodes the Lie generator
struct Advection{T}
	p::T
end

function LinearAlgebra.mul!(f, aD::Advection, U, α, β)
	rmul!(f, β)
	applyD_add!(f, U, aD.p, α)
end

uold = vcat(u0, (1 .- u0))

# we create a TW problem
probTW = BK.TWProblem(prob, Advection(par_cat), copy(uold))

We now define the $U_0$ profile

uold = vcat(u0, (1 .- u0))
Duold = zero(uold); applyD_add!(Duold, uold, par_cat,1)

# update problem parameters for front problem
par_cat_wave = (par_cat..., u0Du0 = dot(uold, Duold), Du0 = Duold, uold = uold)

Let us find the front using newton

front = newton(probTW, vcat(U0, 0.1), NewtonPar(verbose = true))
println("norm front = ", front.u[1:end-1] |> norminf, ", speed = ", front.u[end])

┌─────────────────────────────────────────────────────┐
│ Newton step         residual      linear iterations │
├─────────────┬──────────────────────┬────────────────┤
│       0     │       2.7909e+00     │        0       │
│       1     │       4.1744e-01     │        1       │
│       2     │       6.7151e-02     │        1       │
│       3     │       9.5457e-03     │        1       │
│       4     │       4.0213e-05     │        1       │
│       5     │       8.7134e-10     │        1       │
│       6     │       1.2310e-14     │        1       │
└─────────────┴──────────────────────┴────────────────┘
norm front = 1.0, speed = 0.26271255049274705

plotsol(front.u[1:end-1], title="front solution")

Continuation of front solutions

Following ^[Malham], the modulated fronts are solutions of the following DAE

\[\begin{array}{l}\tag{DAE} u_{t}=a u_{x x}+s\cdot u_x-u f(v)\\ v_{t}=v_{x x}+s\cdot v_x+u f(v)\\ 0 = \left\langle U, \partial_\xi U_0 \right\rangle \end{array}\]

which can be written with a PDE $M_aU_t = G(u)$ with mass matrix $M_a = (Id, Id, 0)$. We have already written the vector field of (MF) in the problem probTW.

Having found a front $U^f$, we can continue it as function of the parameter $a$ and detect instabilities. The stability of the front is linked to the eigenelements $(\lambda, V)$ solution of the generalized eigenvalue problem:

\[\lambda M_a\cdot V = dG(U^f)\cdot V.\]

This is handled automatically when calling continuation on a TWProblem.

optn = NewtonPar(tol = 1e-8)
opt_cont_br = ContinuationPar(p_min = 0.015, p_max = 0.18, newton_options = optn, ds= -0.001, plot_every_step = 2, detect_bifurcation = 3, nev = 10, n_inversion = 6)
br_TW = continuation(probTW, front.u, PALC(), opt_cont_br)
plot(br_TW, legend = :topright)

We have detected a Hopf instability in front dynamics, this will give rise of modulated fronts.

References