improve examples

bifurcationkit · May 8, 2023 · 6f70f4f · 6f70f4f
1 parent 7a7fb3c
commit 6f70f4f
Show file tree

Hide file tree

Showing 4 changed files with 36 additions and 10 deletions.
diff --git a/examples/ikeda.jl b/examples/ikeda.jl
@@ -65,6 +65,13 @@ probpo = PeriodicOrbitTrapProblem(M = 100, jacobian = :DenseAD, N = 1)
 	args_po...,
 	# ampfactor = 1/0.467829783456199 ,#* 0.014,
 	δ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")
+		@show state.x[end]
+		state.it < 18
+	end
 	)
 
 plot(br, br_pocoll)

diff --git a/examples/neuron.jl b/examples/neuron.jl
@@ -65,29 +65,40 @@ opts_po_cont = ContinuationPar(dsmax = 0.1, ds= 0.0001, dsmin = 1e-4, pMax = 10.
 
 	# arguments for periodic orbits
 	args_po = (	recordFromSolution = (x, p) -> begin
-			xtt = BK.getPeriodicOrbit(p.prob, x, nothing)
+			_par = BK.getParams(p.prob)
+			xtt = BK.getPeriodicOrbit(p.prob, x, @set _par.a21=p.p)
 			return (max = maximum(xtt[1,:]),
 					min = minimum(xtt[1,:]),
-					period = getPeriod(p.prob, x, nothing))
+					period = getPeriod(p.prob, x, @set _par.a21=p.p))
 		end,
 		plotSolution = (x, p; k...) -> begin
-			xtt = BK.getPeriodicOrbit(p.prob, x, nothing)
+			_par = BK.getParams(p.prob)
+			xtt = BK.getPeriodicOrbit(p.prob, x, @set _par.a21=p.p)
 			plot!(xtt.t, xtt[1,:]; label = "V1", k...)
 			plot!(xtt.t, xtt[2,:]; label = "V2", k...)
 			plot!(br2; subplot = 1, putspecialptlegend = false)
 			end,
 		normC = norminf)
 
-probpo = PeriodicOrbitOCollProblem(60, 4; N = 2)
-	# probpo = PeriodicOrbitTrapProblem(M = 2000, jacobian = :DenseAD, N = 2)
+probpo = PeriodicOrbitOCollProblem(50, 3; N = 2)
+	# probpo = PeriodicOrbitTrapProblem(M = 150, N = 2; jacobian = :DenseAD)
 	br_pocoll = @time continuation(
 		br2, 1, opts_po_cont,
 		probpo;
+		alg = PALC(tangent = Bordered()),
 		verbosity = 2,	plot = true,
 		args_po...,
-		# ampfactor = 1/0.3593 * 0.0610*2.2,
+		# eigsolver = BK.FloquetCollGEV(DefaultEig(), length(probpo), probpo.N),
 		δp = 0.001,
 		normC = norminf,
+		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")
+			# @show state.x[end]
+			# @show state.f[end]
+			state.it < 6
+		end
 		)
 
 plot(br2, br_pocoll)
@@ -96,7 +107,7 @@ plot(br_pocoll, vars = (:param, :period))
 # plot the periodic orbit
 plot(layout = 2)
 	for ii = 1:10:110
-		solpo = BK.getPeriodicOrbit(br_pocoll.γ.prob.prob, br_pocoll.sol[ii].x, nothing)
+		solpo = BK.getPeriodicOrbit(br_pocoll.γ.prob.prob, br_pocoll.sol[ii].x, 1)
 		plot!(solpo.t ./ solpo.t[end], solpo.u[1,:], label = "", subplot = 1)
 	end
 	xlabel!("t / period", subplot = 1)

diff --git a/examples/neuronV2.jl b/examples/neuronV2.jl
@@ -96,7 +96,7 @@ opts_po_cont = ContinuationPar(dsmax = 0.1, ds= -0.0001, dsmin = 1e-4, pMax = 10
 			end,
 		normC = norminf)
 
-probpo = PeriodicOrbitOCollProblem(140, 3; N = 2)
+probpo = PeriodicOrbitOCollProblem(100, 3; N = 2)
 	# probpo = PeriodicOrbitTrapProblem(M = 2000, jacobian = :DenseAD, N = 2)
 	br_pocoll = @time continuation(
 		br, 1, opts_po_cont,
@@ -108,6 +108,14 @@ probpo = PeriodicOrbitOCollProblem(140, 3; N = 2)
 		ampfactor = 1.42,
 		δp = 0.001,
 		normC = norminf,
+		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")
+			# @show state.x[end]
+			# @show state.f[end]
+			state.it < 16
+		end
 		)
 ################################################################################
 using  DifferentialEquations

diff --git a/src/NormalForms.jl b/src/NormalForms.jl
@@ -81,8 +81,8 @@ end
 
 function BK.hopfNormalForm(prob::SDDDEBifProblem, pt::BK.Hopf, ls; verbose::Bool = false)
 	@error "Normal form for SD-DDE is not implemented"
-	a = Complex{eltype(prob)}(1, 0)
-	b = Complex{eltype(prob)}(1, 0)
+	a = Complex{eltype(pt.x0)}(1, 0)
+	b = Complex{eltype(pt.x0)}(1, 0)
 	pt.nf = (a = a, b = b)
 	if real(a) * real(b) < 0
 		pt.type = :SuperCritical