Pages = ["cgl1dwave.md"]
Depth = 3
We look at the Ginzburg-Landau equations in 1d with periodic boundary condition. The goal of this tutorial is to show how one can study a Hopf bifurcation with symmetry group
The equations are as follows
with periodic boundary conditions. We discretize the circle
using Revise
using DiffEqOperators, ForwardDiff
using BifurcationKit, LinearAlgebra, Plots, SparseArrays, Parameters, Setfield
const BK = BifurcationKit
const FD = ForwardDiff
# plotting utilities
plotsol!(x, m, n; np = n, k...) = heatmap!(reshape(x[1:end-1],m,n)[1:np,:]; color = :viridis, k...)
contoursol!(x, m, n; np = n, k...) = contour!(reshape(x[1:end-1],m,n)[1:np,:]; color = :viridis, k...)
plotsol(x,m,n;k...) = (plot();plotsol!(x,m,n;k...))
function Laplacian1D(Nx, lx)
hx = 2lx/Nx
T = typeof(hx)
D2x = CenteredDifference(2, 2, hx, Nx)
D1x = CenteredDifference(1, 2, hx, Nx)
Qx = PeriodicBC(T)
Δ = sparse(D2x * Qx)[1] |> sparse
D = sparse(D1x * Qx)[1] |> sparse
return Δ, D
end
# add the nonlinearity to f
@views function NL!(f, u, p)
@unpack r, μ, ν, c3, c5 = p
n = div(length(u), 2)
u1 = u[1:n]
u2 = u[n+1:2n]
ua = u1.^2 .+ u2.^2
f1 = f[1:n]
f2 = f[n+1:2n]
f1 .+= @. r * u1 - ν * u2 - ua * (c3 * u1 - μ * u2) - c5 * ua^2 * u1
f2 .+= @. r * u2 + ν * u1 - ua * (c3 * u2 + μ * u1) - c5 * ua^2 * u2
return f
end
# full functional
function Fcgl!(f, u, p, t = 0)
mul!(f, p.Δ, u)
NL!(f, u, p)
end
Fcgl(u, p, t = 0) = Fcgl!(similar(u), u, p)
# analytical expression off the jacobian
@views function Jcgl(u, p)
@unpack r, μ, ν, c3, c5, Δ = p
n = div(length(u), 2)
u1 = u[1:n]
u2 = u[n+1:2n]
ua = u1.^2 .+ u2.^2
f1u = zero(u1)
f2u = zero(u1)
f1v = zero(u1)
f2v = zero(u1)
@. f1u = r - 2 * u1 * (c3 * u1 - μ * u2) - c3 * ua - 4 * c5 * ua * u1^2 - c5 * ua^2
@. f1v = -ν - 2 * u2 * (c3 * u1 - μ * u2) + μ * ua - 4 * c5 * ua * u1 * u2
@. f2u = ν - 2 * u1 * (c3 * u2 + μ * u1) - μ * ua - 4 * c5 * ua * u1 * u2
@. f2v = r - 2 * u2 * (c3 * u2 + μ * u1) - c3 * ua - 4 * c5 * ua * u2 ^2 - c5 * ua^2
jacdiag = vcat(f1u, f2v)
Δ + spdiagm(0 => jacdiag, n => f1v, -n => f2u)
end
nothing #hide
We then define a problem for computing the bifurcations of the trivial state
# space discretization
n = 50
l = pi
Δ, D = Laplacian1D(n, l)
# model parameters
par_cgl = (r = 0.0, μ = 0.5, ν = 1.0, c3 = -1.0, c5 = 1.0, Δ = blockdiag(Δ, Δ), Db = blockdiag(D, D), δ = 1.0, N = 2n)
# initial guess
sol0 = zeros(par_cgl.N)
# bifurcation problem
prob = BifurcationProblem(Fcgl, sol0, par_cgl, (@lens _.r); J = Jcgl)
opt_newton = NewtonPar(tol = 1e-9, max_iterations = 20)
opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.15, ds = 0.001, p_max = 2.5,
detect_bifurcation = 3, nev = 9, newton_options = (@set opt_newton.verbose = false), max_steps = 100, n_inversion = 8, max_bisection_steps = 20)
br = continuation(prob, PALC(), opts_br, verbosity = 0)
The first bifurcation point is a regular Hopf bifurcation in the zero mode, i.e.
We focus on the br
), with frequency BifurcationKit.jl
. We write
By the center manifold theory1, one has
Using the normal form, one finds standing waves
function guessFromHopfO2(branch, ind_hopf, eigsolver, M, A, B = 0.; phase = 0, k = 1.)
specialpoint = branch.specialpoint[ind_hopf]
# parameter value at the Hopf point
p_hopf = specialpoint.param
# frequency at the Hopf point
ωH = imag(branch.eig[specialpoint.idx].eigenvals[specialpoint.ind_ev]) |> abs
# eigenvectors for the eigenvalues iω
ζ0 = geteigenvector(eigsolver, br.eig[specialpoint.idx][2], specialpoint.ind_ev)
ζ0 ./= norm(ζ0)
ζ1 = geteigenvector(eigsolver, br.eig[specialpoint.idx][2], specialpoint.ind_ev - 2)
ζ1 ./= norm(ζ1)
orbitguess = [real.(specialpoint.x .+
A .* ζ0 .* exp(2pi * complex(0, 1) .* (ii/(M-1) - phase)) .+
B .* ζ1 .* exp(2pi * complex(0, 1) .* (ii/(M-1) - phase))) for ii in 0:M-1]
return (; p = p_hopf, period = 2pi/ωH, guess = orbitguess, x0 = specialpoint.x, ζ0 = ζ0, ζ1 = ζ1)
end
nothing #hide
We can use this function to effectively build a guess for the travelling wave:
M = 50 # number of time slices (plotting purposes)
r_hopf, Th, orbitguess2, hopfpt, eigvec = guessFromHopfO2(br, 2, opt_newton.eigsolver, M, 1. + 0.0im, 1+0.0im; k = 2.) #TW
uold = copy(orbitguess2[1][1:2n])
# we create a TW problem
probTW = TWProblem(re_make(prob, params = setproperties(par_cgl; r = r_hopf - 0.01)), par_cgl.Db, uold; jacobian = :FullLU)
# refine the guesss
wave = newton(probTW, vcat(uold, 0),
NewtonPar(verbose = true, max_iterations = 50),
)
println("norm wave = ", wave.u[1:end-1] |> norminf)
plot(wave.u[1:end-1]; linewidth = 5, label = "solution")
plot!(uold, color = :blue, label="guess")
Note that in the following code, a generalized eigensolver is automatically created during the call to continuation
which properly computes the stability of the wave.
amplitude(x) = maximum(x) - minimum(x)
optn = NewtonPar(tol = 1e-8, verbose = true, max_iterations = 10)
opt_cont_br = ContinuationPar(p_min = 0.015, p_max = 2.5, newton_options = optn, ds= 0.001, dsmax = 0.1, detect_bifurcation = 3, nev = 10, max_steps = 190, n_inversion = 6)
br_TW = @time continuation(probTW, wave.u, PALC(), opt_cont_br;
record_from_solution = (x, p) -> (u∞ = maximum(x[1:n]), s = x[end], amp = amplitude(x[1:n])),
plot_solution = (x, p; k...) -> (plot!(x[1:end-1];k...);plot!(br,subplot=1, legend=false)),
finalise_solution = (z, tau, step, contResult; k...) -> begin
amplitude(z.u[n+1:2n]) > 0.01
end, bothside = true)
plot(br, br_TW, legend = :bottomright, branchlabel =["","TW"])
We note that the branch of travelling wave solutions has a Hopf bifurcation point at which point Modulated Travelling waves will emerge. This will be analyzed in the future.
Footnotes
-
Haragus, Mariana, and Gérard Iooss. Local Bifurcations, Center Manifolds, and Normal Forms in Infinite-Dimensional Dynamical Systems. London: Springer London, 2011. https://doi.org/10.1007/978-0-85729-112-7.