using RecursiveArrayTools # for bifurcation point handling
import Base: show
###############################################################################################
# Structure to hold result
@with_kw struct ContResult{T, vectype, eigenvectype}
# this vector is used to hold (param, printsolution(u), Newton iterations, ds)
branch::VectorOfArray{T, 2, Array{Vector{T}, 1}} = VectorOfArray([zeros(T, 4)])
# the following holds the eigen elements at the index of the point along the curve. This index is the last element of eig[1] (for example). Recording this index is useful for storing only some eigenelements and not all of them along the curve
eig::Vector{Tuple{Array{Complex{T}, 1}, eigenvectype, Int64}}
# the following holds information about the detected bifurcation points like
# [(:none, step, printsolution(u), u, tau, eigenvalue_index)] where tau is the tangent along the curve and eigenvalue_index is the index of the eigenvalue in eig (see above) which change stability
bifpoint::Vector{Tuple{Symbol, Int64, T, T, vectype, vectype, Int64}}
# whether the associated point is linearly stable
stability::Vector{Bool} = [false]
# number of eigenvalues with positive real part and non zero imaginary part
n_imag::Vector{Int64}
# number of eigenvalues with positive real part
n_unstable::Vector{Int64}
end
function show(io::IO, br::PseudoArcLengthContinuation.ContResult)
println(io, "Branch number of points: ", length(br.branch))
if length(br.bifpoint) >0
println(io, "Bifurcation points:")
for ii in eachindex(br.bifpoint)
bp = br.bifpoint[ii]
println(io, "- $ii, ", bp[1], " point, at p = ", bp[3])
end
end
end
###############################################################################################
function displayIteration(i, funceval, residual, itlinear = 0)
(i == 0) && println("\n Newton Iterations \n Iterations Func-count f(x) Linear-Iterations\n")
if length(itlinear)==1
@printf("%9d %16d %14.4e %9d\n", i, funceval, residual, itlinear);
else
@printf("%9d %16d %14.4e (%6d, %6d)\n", i, funceval, residual, itlinear[1], itlinear[2]);
end
end
###############################################################################################
"""
Plot the continued branch of solutions
"""
function plotBranchCont(contres::ContResult{T}, sol::M, contparms, plotuserfunction) where {T, vectype, M<:BorderedVector{vectype, T}}
colorbif = Dict(:fold => :black, :hopf => :red, :bp => :blue, :nd => :magenta)
branch = contres.branch
if contparms.computeEigenValues == false
l = Plots.@layout [a{0.5w} b{0.5w}; c d]
else
l = Plots.@layout [a{0.5w} b{0.5w}; c d;e]
end
Plots.plot(layout=l)
# plot the branch of solutions
plotBranch!(contres, xlabel="p", ylabel="|x|", label="", subplot=1 )
plot!(branch[1, :], xlabel="s", ylabel="p", label="", subplot=2)
plot!(branch[2, :], xlabel="it", ylabel="|x|", label="", subplot=3)
# add the bifurcation points along the branch
if length(contres.bifpoint)>1
scatter!(map(x->x[2],contres.bifpoint[2:end]), map(x->x[4],contres.bifpoint[2:end]), label="", color = map(x->colorbif[x[1]],contres.bifpoint[2:end]), markersize=3, markerstrokewidth=0, subplot=3) |> display
end
if contparms.computeEigenValues
for ii=1:length(contres.eig)
scatter!(real.(contres.eig[ii][1]), imag.(contres.eig[ii][1]), subplot=5, label="", markersize = 1, color=:black)
end
end
plotuserfunction(sol.u)
display(title!(""))
end
"""
Plot the branch of solutions from a `ContResult`. You can also pass parameters like `plotBranch(br, marker = :dot)`.
For the continuation diagram, the legend is as follows `(:fold => :black, :hopf => :red, :bp => :blue, :nd => :magenta, :none => :yellow)`
"""
function plotBranch(contres::ContResult; kwargs...)
plot([],[],label="")
plotBranch!(contres; kwargs...)
end
"""
Plot all the branches contained in `brs` in a single figure. Convenient when many bifurcation diagram have been computed.
"""
function plotBranch(brs::Vector; kwargs...)
plotBranch(brs[1]; kwargs...)
for ii=2:length(brs)
plotBranch!(brs[ii]; kwargs...) |> display
end
end
"""
Append to the current plot the plot of the branch of solutions from a `ContResult`. You can also pass parameters like `plotBranch!(br, marker = :dot)`
"""
function plotBranch!(contres::ContResult; kwargs...)
colorbif = Dict(:fold => :black, :hopf => :red, :bp => :blue, :nd => :magenta, :none => :yellow)
branch = contres.branch
if length(contres.stability) > 2
plot!(branch[1, :], branch[2, :], linewidth = 1 .+ 3contres.stability ; kwargs...)
else
plot!(branch[1, :], branch[2, :]; kwargs...)
end
# add the bifurcation points along the branch
if length(contres.bifpoint)>=1
id = 1
contres.bifpoint[1][1] == :none ? id = 2 : id = 1
scatter!(map(x->x[3],contres.bifpoint[id:end]), map(x->x[4],contres.bifpoint[id:end]), label="", color = map(x->colorbif[x[1]],contres.bifpoint[id:end]), markersize=3, markerstrokewidth=0 ; kwargs...) |> display
end
end
###############################################################################################
function detectBifucation(contparams, contResult, z, tau, printsolution, verbosity)
branch = contResult.branch
# Fold point detection based on continuation parameter monotony
if contparams.detect_fold && size(branch)[2] > 2 && (branch[1, end] - branch[1, end-1])* (branch[1, end-1] - branch[1, end-2]) < 0
(verbosity > 1) && printstyled(color=:red, "Fold bifurcation point!! between $(branch[1, end-1]) and $(branch[1, end]) \n")
push!(contResult.bifpoint, (:fold,
length(branch)-1,
branch[1, end-1],
printsolution(z.u), copy(z.u), copy(tau.u) / norm(tau.u), 0))
end
# update number of unstable eigenvalues
n_unstable = mapreduce(x -> round(real(x), digits=6) > 0, +, contResult.eig[end][1])
push!(contResult.n_unstable, n_unstable)
# update number of unstable eigenvalues with nonzero imaginary part
n_imag = mapreduce(x -> (abs(round(imag(x), digits = 6)) > 0) * (round(real(x), digits = 6) > 0), +, contResult.eig[end][1])
push!(contResult.n_imag, n_imag)
# computation of the index of the bifurcating eigenvalue
ind_bif = n_unstable
if n_unstable < contResult.n_unstable[end-1]
ind_bif += 1
end
# Hopf / BP bifurcation point detection based on eigenvalue distribution
if size(branch)[2] > 1
if abs(contResult.n_unstable[end] - contResult.n_unstable[end-1]) == 1
push!(contResult.bifpoint, (:bp,
length(branch)-1,
branch[1, end-1],
printsolution(z.u),
z.u,
tau.u ./ norm(tau.u), ind_bif))
elseif abs(contResult.n_unstable[end] - contResult.n_unstable[end-1]) == 2
if abs(contResult.n_imag[end] - contResult.n_imag[end-1]) == 2
push!(contResult.bifpoint, (:hopf,
length(branch)-1,
branch[1, end-1],
printsolution(z.u),
copy(z.u), 0tau.u, ind_bif))
else
push!(contResult.bifpoint, (:bp,
length(branch)-1,
branch[1, end-1],
printsolution(z.u),
copy(z.u),
copy(tau.u) / norm(tau.u), n_unstable))
end
elseif abs(contResult.n_unstable[end] - contResult.n_unstable[end-1]) >2
push!(contResult.bifpoint, (:nd,
length(branch)-1,
branch[1, end-1],
printsolution(z.u),
copy(z.u),
copy(tau.u) / norm(tau.u), ind_bif))
end
end
end
###############################################################################################
"""
Compute Jacobian by Finite Differences
"""
function finiteDifferences(F, x::AbstractVector)
f = F(x)
epsi = 1e-9
N = length(x)
J = zeros(N, N)
x1 = copy(x)
for i=1:N
x1[i] += epsi
J[:, i] .= (F(x1) .- F(x)) / epsi
x1[i] -= epsi
end
return J
end
###############################################################################################
"""
Save solution / data in JLD2 file
- `filename` is for example "example.jld2"
- `sol` is the solution
- `p` is the parameter
- `i` is the index of the solution to be saved
"""
function saveSolution(filename, sol, p, i::Int64, br::ContResult, contParam)
# create a group in the JLD format
jldopen(filename*".jld2", "a+") do file
mygroup = JLD2.Group(file, "sol-$i")
mygroup["sol"] = sol
mygroup["param"] = p
end
jldopen(filename*"-branch.jld2", "w") do file
file["branch"] = br
file["contParam"] = contParam
end
end
###############################################################################################
# this trick is extracted from KrylovKit. It allows for the Jacobian to be specified as a matrix (sparse / dense) or as a function.
apply(A::AbstractMatrix, x::AbstractVector) = A * x
apply(f, x) = f(x)
apply!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, A, x)
apply!(y, f, x) = copyto!(y, f(x))