test periodic orbits for sdde

remove unused code in Hopf.jl

src/Problems.jl

@@ -51,6 +51,7 @@ BK.isSymmetric(::ConstantDDEBifProblem) = false
 BK.getVectorType(prob::ConstantDDEBifProblem{Tvf, Tdf, Tu, Td, Tp, Tl, Tplot, Trec}) where {Tvf, Tdf, Tu, Td, Tp, Tl <: Lens, Tplot, Trec} = Tu
 BK.getLens(prob::ConstantDDEBifProblem) = prob.lens
 BK.hasAdjoint(prob::ConstantDDEBifProblem) = true
+BK.hasAdjointMF(prob::ConstantDDEBifProblem) = false
 BK.getDelta(prob::ConstantDDEBifProblem) = prob.δ
 BK.d2F(prob::ConstantDDEBifProblem, x, p, dx1, dx2) = BK.d2F(prob.VF, x, p, dx1, dx2)
 BK.d3F(prob::ConstantDDEBifProblem, x, p, dx1, dx2, dx3) = BK.d3F(prob.VF, x, p, dx1, dx2, dx3)

src/codim2/Hopf.jl

@@ -320,221 +320,6 @@ function BK.continuationHopf(prob::AbstractDDEBifurcationProblem,
 					kwargs...)
 end
 
-function BK.continuationHopf(prob_vf::AbstractDDEBifurcationProblem, alg::BK.AbstractContinuationAlgorithm,
-				hopfpointguess::ArrayPartition, par,
-				lens1::Lens, lens2::Lens,
-				eigenvec, eigenvec_ad,
-				options_cont::ContinuationPar ;
-				updateMinAugEveryStep = 0,
-				normC = norm,
-				bdlinsolver::BK.AbstractBorderedLinearSolver = MatrixBLS(),
-				jacobian_ma::Symbol = :autodiff,
-				computeEigenElements = false,
-				usehessian = true,
-				massmatrix = LinearAlgebra.I,
-				kwargs...) where {Tb, vectype}
-	@assert lens1 != lens2 "Please choose 2 different parameters. You only passed $lens1"
-	@assert lens1 == BK.getLens(prob_vf)
-
-	# options for the Newton Solver inheritated from the ones the user provided
-	options_newton = options_cont.newtonOptions
-	threshBT = 100options_newton.tol
-
-	hopfPb = HopfDDEProblem(
-		prob_vf,
-		BK._copy(eigenvec_ad),	# this is a ≈ null space of (J - iω I)^*
-		BK._copy(eigenvec), 	# this is b ≈ null space of  J - iω I
-		options_newton.linsolver,
-		# do not change linear solver if user provides it
-		@set bdlinsolver.solver = (isnothing(bdlinsolver.solver) ? options_newton.linsolver : bdlinsolver.solver);
-		usehessian = usehessian,
-		massmatrix = massmatrix)
-
-	# Jacobian for the Hopf problem
-	if jacobian_ma == :autodiff
-		# hopfpointguess = vcat(hopfpointguess.u, hopfpointguess.p)
-		prob_h = BK.HopfMAProblem(hopfPb, BK.AutoDiff(), hopfpointguess, par, lens2, prob_vf.plotSolution, prob_vf.recordFromSolution)
-		opt_hopf_cont = @set options_cont.newtonOptions.linsolver = DefaultLS()
-	elseif jacobian_ma == :finiteDifferencesMF
-		hopfpointguess = vcat(hopfpointguess.u, hopfpointguess.p)
-		prob_h = FoldMAProblem(hopfPb, FiniteDifferences(), hopfpointguess, par, lens2, prob_vf.plotSolution, prob_vf.recordFromSolution)
-		opt_hopf_cont = @set options_cont.newtonOptions.linsolver = options_cont.newtonOptions.linsolver
-	else
-		prob_h = HopfMAProblem(hopfPb, nothing, hopfpointguess, par, lens2, prob_vf.plotSolution, prob_vf.recordFromSolution)
-		opt_hopf_cont = @set options_cont.newtonOptions.linsolver = HopfLinearSolverMinAug()
-	end
-
-	# this functions allows to tackle the case where the two parameters have the same name
-	lenses = BK.getLensSymbol(lens1, lens2)
-
-	# current lyapunov coefficient
-	eTb = eltype(hopfpointguess)
-	hopfPb.l1 = Complex{eTb}(0, 0)
-	hopfPb.BT = one(eTb)
-	hopfPb.GH = one(eTb)
-
-	# this function is used as a Finalizer
-	# it is called to update the Minimally Augmented problem
-	# by updating the vectors a, b
-	function updateMinAugHopf(z, tau, step, contResult; kUP...)
-		# we first check that the continuation step was successful
-		# if not, we do not update the problem with bad information!
-		success = get(kUP, :state, nothing).converged
-		(~modCounter(step, updateMinAugEveryStep) || success == false) && return true
-		x = getVec(z.u, hopfPb)	# hopf point
-		p1, ω = getP(z.u, hopfPb)
-		p2 = z.p		# second parameter
-		newpar = set(par, lens1, p1)
-		newpar = set(newpar, lens2, p2)
-
-		a = hopfPb.a
-		b = hopfPb.b
-
-		# expression of the jacobian
-		J_at_xp = jacobian(hopfPb.prob_vf, x, newpar)
-
-		# compute new b
-		T = typeof(p1)
-		newb = hopfPb.linbdsolver(J_at_xp, a, b, T(0), hopfPb.zero, T(1); shift = Complex(0, -ω), Mass = hopfPb.massmatrix)[1]
-
-		# compute new a
-		JAd_at_xp = hasAdjoint(hopfPb) ? jad(hopfPb.prob_vf, x, newpar) : adjoint(J_at_xp)
-		newa = hopfPb.linbdsolver(JAd_at_xp, b, a, T(0), hopfPb.zero, T(1); shift = Complex(0, ω), Mass = adjoint(hopfPb.massmatrix))[1]
-
-		hopfPb.a .= newa ./ normC(newa)
-		# do not normalize with dot(newb, hopfPb.a), it prevents BT detection
-		hopfPb.b .= newb ./ normC(newb)
-
-		# we stop continuation at Bogdanov-Takens points
-
-		# CA NE DEVRAIT PAS ETRE ISSNOT?
-		isbt = isnothing(contResult) ? true : isnothing(findfirst(x -> x.type in (:bt, :ghbt, :btgh), contResult.specialpoint))
-
-		# if the frequency is null, this is not a Hopf point, we halt the process
-		if abs(ω) < threshBT
-			@warn "[Codim 2 Hopf - Finalizer] The Hopf curve seems to be close to a BT point: ω ≈ $ω. Stopping computations at ($p1, $p2). If the BT point is not detected, try lowering Newton tolerance or dsmax."
-		end
-
-		# call the user-passed finalizer
-		finaliseUser = get(kwargs, :finaliseSolution, nothing)
-		resFinal = isnothing(finaliseUser) ? true : finaliseUser(z, tau, step, contResult; prob = hopfPb, kUP...)
-
-		return abs(ω) >= threshBT && isbt && resFinal
-	end
-
-	function testBT_GH(iter, state)
-		z = getx(state)
-		x = getVec(z, hopfPb)		# hopf point
-		p1, ω = getP(z, hopfPb)		# first parameter
-		p2 = getp(state)			# second parameter
-		newpar = set(par, lens1, p1)
-		newpar = set(newpar, lens2, p2)
-
-		probhopf = iter.prob.prob
-
-		a = probhopf.a
-		b = probhopf.b
-
-		# expression of the jacobian
-		J_at_xp = BK.jacobian(probhopf.prob_vf, x, newpar)
-
-		# compute eigenvector
-		T = typeof(p1)
-		ζ = @. z.x[2] + im * z.x[3]
-		ζ ./= normC(ζ)
-
-		# test function for Bogdanov-Takens
-		probhopf.BT = ω
-
-		return probhopf.BT, probhopf.GH
-	end
-
-	# the following allows to append information specific to the codim 2 continuation to the user data
-	_printsol = get(kwargs, :recordFromSolution, nothing)
-	_printsol2 = isnothing(_printsol) ?
-		(u, p; kw...) -> (; zip(lenses, (getP(u, hopfPb)[1], p))..., ω = getP(u, hopfPb)[2], l1 = hopfPb.l1, BT = hopfPb.BT, GH = hopfPb.GH, BK.namedprintsol(BK.recordFromSolution(prob_vf)(getVec(u, hopfPb), p; kw...))...) :
-		(u, p; kw...) -> (; BK.namedprintsol(_printsol(getVec(u, hopfPb), p; kw...))..., zip(lenses, (getP(u, hopfPb)[1], p))..., ω = getP(u, hopfPb)[2], l1 = hopfPb.l1, BT = hopfPb.BT, GH = hopfPb.GH)
-
-	prob_h = reMake(prob_h, recordFromSolution = _printsol2)
-
-	# eigen solver
-	eigsolver = HopfDDEEig(BK.getsolver(opt_hopf_cont.newtonOptions.eigsolver))
-
-	# event for detecting codim 2 points
-	if computeEigenElements
-		event = PairOfEvents(ContinuousEvent(2, testBT_GH, computeEigenElements, ("bt", "gh"), threshBT), BifDetectEvent)
-		# careful here, we need to adjust the tolerance for stability to avoid
-		# spurious ZH or HH bifurcations
-		@set! opt_hopf_cont.tolStability = 10opt_hopf_cont.newtonOptions.tol
-	else
-		event = ContinuousEvent(2, testBT_GH, false, ("bt", "gh"), threshBT)
-	end
-
-	prob_h = reMake(prob_h, recordFromSolution = _printsol2)
-
-	# solve the hopf equations
-	br = continuation(
-		prob_h, alg,
-		(@set opt_hopf_cont.newtonOptions.eigsolver = eigsolver);
-		kwargs...,
-		kind = BK.HopfCont(),
-		linearAlgo = BorderingBLS(solver = opt_hopf_cont.newtonOptions.linsolver, checkPrecision = false),
-		normC = normC,
-		finaliseSolution = updateMinAugEveryStep ==0 ? get(kwargs, :finaliseSolution, BK.finaliseDefault) : updateMinAugHopf,
-		event = event
-	)
-	@assert ~isnothing(br) "Empty branch!"
-	return correctBifurcation(br)
-end
-
-# function BK.continuationHopf(prob::AbstractDDEBifurcationProblem,
-# 						br::BK.AbstractBranchResult, ind_hopf::Int64,
-# 						lens2::Lens,
-# 						options_cont::ContinuationPar = br.contparams;
-# 						alg = br.alg,
-# 						startWithEigen = false,
-# 						normC = norm,
-# 						kwargs...)
-# 	hopfpointguess = HopfPoint(br, ind_hopf)
-# 	ω = hopfpointguess.p[2]
-# 	bifpt = br.specialpoint[ind_hopf]
-
-# 	@assert ~isnothing(br.eig) "The branch contains no eigen elements. This is strange because a Hopf point was detected. Please open an issue on the website."
-
-# 	@assert ~isnothing(br.eig[1].eigenvecs) "The branch contains no eigenvectors for the Hopf point. Please provide one."
-
-# 	ζ = geteigenvector(br.contparams.newtonOptions.eigsolver, br.eig[bifpt.idx].eigenvecs, bifpt.ind_ev)
-# 	ζ ./= normC(ζ)
-# 	ζad = conj.(ζ)
-
-# 	p = bifpt.param
-# 	parbif = BK.setParam(br, p)
-
-# 	hopfpointguess = ArrayPartition(hopfpointguess.u, real(ζ), imag(ζ), hopfpointguess.p)
-
-# 	if startWithEigen
-# 		# computation of adjoint eigenvalue
-# 		λ = Complex(0, ω)
-# 		# jacobian at bifurcation point
-# 		L = BK.jacobian(prob, bifpt.x, parbif)
-
-# 		# jacobian adjoint at bifurcation point
-# 		_Jt = ~BK.hasAdjoint(prob) ? adjoint(L) : jad(prob, bifpt.x, parbif)
-
-# 		ζstar, λstar = BK.getAdjointBasis(_Jt, conj(λ), br.contparams.newtonOptions.eigsolver; nev = br.contparams.nev, verbose = false)
-# 		ζad .= ζstar ./ dot(ζstar, ζ)
-# 	end
-
-# 	return BK.continuationHopf(br.prob, alg,
-# 					hopfpointguess, parbif,
-# 					BK.getLens(br), lens2,
-# 					ζ, ζad,
-# 					options_cont ;
-# 					normC = normC,
-# 					kwargs...)
-# end
-
-
 # structure to compute the eigenvalues along the Hopf branch
 struct HopfDDEEig{S} <: BK.AbstractCodim2EigenSolver
 	eigsolver::S

src/codim2/codim2.jl

@@ -44,10 +44,10 @@ for op in (:HopfDDEProblem,)
 			massmatrix::Tmass
 		end
 
-		@inline hasHessian(pb::$op) = hasHessian(pb.prob_vf)
+		@inline BK.hasHessian(pb::$op) = BK.hasHessian(pb.prob_vf)
 		@inline BK.isSymmetric(pb::$op) = false
-		@inline hasAdjoint(pb::$op) = hasAdjoint(pb.prob_vf)
-		@inline hasAdjointMF(pb::$op) = hasAdjointMF(pb.prob_vf)
+		@inline BK.hasAdjoint(pb::$op) = BK.hasAdjoint(pb.prob_vf)
+		@inline BK.hasAdjointMF(pb::$op) = BK.hasAdjointMF(pb.prob_vf)
 		@inline BK.isInplace(pb::$op) = BK.isInplace(pb.prob_vf)
 		@inline BK.getLens(pb::$op) = BK.getLens(pb.prob_vf)
 		jad(pb::$op, args...) = jad(pb.prob_vf, args...)

test/codim2.jl

@@ -42,6 +42,12 @@ brhopf = continuation(br, 1, (@lens _.c),
 			bothside = true,
 			startWithEigen = true)
 
+BK.isInplace(brhopf.prob)
+BK.isSymmetric(brhopf.prob)
+BK.hasAdjoint(brhopf.prob.prob)
+BK.hasAdjointMF(brhopf.prob.prob)
+BK.getLens(brhopf.prob.prob)
+
 brhopf2 = continuation(br, 2, (@lens _.c),
 			setproperties(br.contparams, detectBifurcation = 1, dsmax = 0.01, maxSteps = 100, pMax = 1.1, pMin = -0.1,ds = -0.01);
 			verbosity = 0, plot = false,

test/normalform.jl

@@ -51,6 +51,8 @@ BK.isInplace(prob)
 BK.isSymmetric(prob)
 BK.getVectorType(prob)
 show(prob)
+
+DDEBifurcationKit.newtonHopf(br, 2)
 ##########################################################################################
 # SDDE, test dummy Hopf normal form
 function humpriesVF(x, xd, p)

test/pocoll-sdde.jl

@@ -0,0 +1,80 @@
+# using Revise, Plots
+using DDEBifurcationKit, Parameters, Setfield, LinearAlgebra
+using BifurcationKit
+const BK = BifurcationKit
+const DDEBK = DDEBifurcationKit
+
+# sup norm
+norminf(x) = norm(x, Inf)
+
+function humpriesVF(x, xd, p)
+   @unpack κ1,κ2,γ,a1,a2,c = p
+   [
+      -γ * x[1] - κ1 * xd[1][1] - κ2 * xd[2][1]
+   ]
+end
+
+function delaysF(x, par)
+   [
+      par.a1 + par.c * x[1],
+      par.a2 + par.c * x[1],
+   ]
+end
+
+
+pars = (κ1=0.,κ2=2.3,a1=1.3,a2=6,γ=4.75,c=1.)
+x0 = zeros(1)
+
+prob = SDDDEBifProblem(humpriesVF, delaysF, x0, pars, (@lens _.κ1))
+
+optn = NewtonPar(verbose = false, eigsolver = DDE_DefaultEig())
+opts = ContinuationPar(pMax = 13., pMin = 0., newtonOptions = optn, ds = -0.01, detectBifurcation = 3, nev = 3, )
+
+alg = PALC()
+br = continuation(prob, alg, opts; verbosity = 0, plot = false, bothside = true)
+
+# plot(br)
+################################################################################
+# computation of periodic orbit
+# continuation parameters
+opts_po_cont = ContinuationPar(dsmax = 0.05, ds= 0.001, dsmin = 1e-4, pMax = 12., pMin=-5., maxSteps = 3,
+	nev = 3, tolStability = 1e-8, detectBifurcation = 0, plotEveryStep = 20, saveSolEveryStep=1)
+	@set! opts_po_cont.newtonOptions.tol = 1e-9
+	@set! opts_po_cont.newtonOptions.verbose = true
+
+	# arguments for periodic orbits
+	args_po = (	recordFromSolution = (x, p) -> begin
+			xtt = BK.getPeriodicOrbit(p.prob, x, nothing)
+			_max = maximum(xtt[1,:])
+			_min = minimum(xtt[1,:])
+			return (amp = _max - _min,
+					max = _max,
+					min = _min,
+					period = getPeriod(p.prob, x, nothing))
+		end,
+		plotSolution = (x, p; k...) -> begin
+			xtt = BK.getPeriodicOrbit(p.prob, x, nothing)
+			plot!(xtt.t, xtt[1,:]; label = "x", k...)
+			plot!(br; subplot = 1, putspecialptlegend = false)
+			end,
+		normC = norminf)
+
+probpo = PeriodicOrbitOCollProblem(200, 2; N = 1)
+	br_pocoll = @time continuation(
+	br, 2, opts_po_cont,
+	probpo;
+	alg = PALC(tangent = Bordered()),
+	# regular continuation options
+	# verbosity = 2,	plot = true,
+	args_po...,
+	ampfactor = 1/0.467829783456199 * 0.1,
+	δp = 0.01,
+	callbackN = (state; k...) -> begin
+		xtt = BK.getPeriodicOrbit(probpo,state.x,nothing)
+		# plot(xtt.t, xtt[1,:], title = "it = $(state.it)") |> display
+		printstyled(color=:red, "amp = ", BK.amplitude(xtt[:,:],1),"\n")
+		printstyled(color=:green, "T = ", (state.x[end]),"\n")
+		@show state.x[end]
+		state.it < 15 && BK.cbMaxNorm(10.0)(state; k...)
+	end
+	)

test/runtests.jl

@@ -5,4 +5,5 @@ using Test
     include("normalform.jl")
     include("codim2.jl")
     include("pocoll.jl")
+    include("pocoll-sdde.jl")
 end