Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Enable setting Δt in minimal_period #25

Merged
merged 5 commits into from
Oct 3, 2024
Merged
Show file tree
Hide file tree
Changes from 3 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 22 additions & 9 deletions src/minimal_period.jl
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,9 @@ period is also called prime, principal or fundamental period.
* `maxiter = 40` : Maximum number of Poincare map iterations. Continuous-time systems only.
If the number of Poincare map iterations exceeds `maxiter`, but the point `u0` has not
returned to `atol` neighborhood of itself, the original period `po.T` is returned.
* `Δt = missing` : The time step between points in the trajectory `minT_po.points`. If `Δt`
is `missing`, then `Δt=minT_po.T/$(default_Δt_partition)` is used. Continuous-time
systems only.

## Description

Expand All @@ -31,25 +34,35 @@ this hyperplane. Using the Poincare map, the hyperplane crossings are checked. T
first crossing that is within `atol` distance of the initial point `u0` is the minimal
period. At most `maxiter` crossings are checked.
"""
function minimal_period(ds::DynamicalSystem, po::PeriodicOrbit; kwargs...)
function minimal_period(
ds::DynamicalSystem,
po::PeriodicOrbit;
Δt=missing,
kwargs...
)
type1 = isdiscretetime(ds)
type2 = isdiscretetime(po)
if type1 == type2
newT = _minimal_period(ds, po.points[1], po.T; kwargs...)
return _set_period(ds, po, newT)
return _set_period(ds, po, newT, Δt)
else
throw(ArgumentError("Both the periodic orbit and the dynamical system have to be either discrete or continuous."))
end
end

function _set_period(ds::DynamicalSystem, po, newT)
if newT == po.T
return po
else
# ensure that continuous po has the same amount of points
isdiscretetime(ds) ? Δt = 1 : Δt = newT/default_Δt_partition
return PeriodicOrbit(complete_orbit(ds, po.points[1], newT; Δt=Δt), newT, po.stable)
function _set_period(ds::DynamicalSystem, po, newT, Δt)
newT == po.T && return po

Dt = 1
if !isdiscretetime(ds)
if ismissing(Δt)
# ensure that continuous po has the same amount of points
Dt = newT/default_Δt_partition
else
Dt = Δt
end
end
return PeriodicOrbit(complete_orbit(ds, po.points[1], newT; Δt=Dt), newT, po.stable)
end

function _minimal_period(ds::DiscreteTimeDynamicalSystem, u0, T; atol=1e-4)
Expand Down
6 changes: 5 additions & 1 deletion test/minimal_period.jl
Original file line number Diff line number Diff line change
Expand Up @@ -31,9 +31,13 @@ end
n = 10
po = PeriodicOrbit(ds, u0, n*T, 0.01)
minT_po = minimal_period(ds, po)
@test length(minT_po.points) == minT_po.T / (minT_po.T / PeriodicOrbits.default_Δt_partition)
@test length(minT_po.points) == PeriodicOrbits.default_Δt_partition
@test (ismissing(po.stable) && ismissing(minT_po.stable)) || (po.stable == minT_po.stable)
@test isapprox(T, minT_po.T; atol=1e-4)

Dt = 0.01
minT_po = minimal_period(ds, po; Δt=Dt)
@test length(minT_po.points) == floor(minT_po.T / Dt)
end

function normalhopf(u, p, t)
Expand Down
Loading