/
poincaremap.jl
287 lines (250 loc) · 10.7 KB
/
poincaremap.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
include("hyperplane.jl")
import Roots
export PoincareMap, current_crossing_time
@deprecate poincaremap PoincareMap
const ROOTS_ALG = Roots.A42()
###########################################################################################
# Type definition
###########################################################################################
"""
PoincareMap <: DiscreteTimeDynamicalSystem
PoincareMap(ds::CoupledODEs, plane; kwargs...) → pmap
A discrete time dynamical system that produces iterations over the Poincaré map[^DatserisParlitz2022]
of the given continuous time `ds`. This map is defined as the sequence of points on the
Poincaré surface of section, which is defined by the `plane` argument.
See also [`StroboscopicMap`](@ref), [`poincaresos`](@ref).
## Keyword arguments
* `direction = -1`: Only crossings with `sign(direction)` are considered to belong to
the surface of section. Negative direction means going from less than ``b``
to greater than ``b``.
* `u0 = nothing`: Specify an initial state.
* `rootkw = (xrtol = 1e-6, atol = 1e-8)`: A `NamedTuple` of keyword arguments
passed to `find_zero` from [Roots.jl](https://github.com/JuliaMath/Roots.jl).
* `Tmax = 1e3`: The argument `Tmax` exists so that the integrator can terminate instead
of being evolved for infinite time, to avoid cases where iteration would continue
forever for ill-defined hyperplanes or for convergence to fixed points,
where the trajectory would never cross again the hyperplane.
If during one `step!` the system has been evolved for more than `Tmax`,
then `step!(pmap)` will terminate and error.
## Description
The Poincaré surface of section is defined as sequential transversal crossings a trajectory
has with any arbitrary manifold, but here the manifold must be a hyperplane.
`PoincareMap` iterates over the crossings of the section.
If the state of `ds` is ``\\mathbf{u} = (u_1, \\ldots, u_D)`` then the
equation defining a hyperplane is
```math
a_1u_1 + \\dots + a_Du_D = \\mathbf{a}\\cdot\\mathbf{u}=b
```
where ``\\mathbf{a}, b`` are the parameters of the hyperplane.
In code, `plane` can be either:
* A `Tuple{Int, <: Real}`, like `(j, r)`: the plane is defined
as when the `j`th variable of the system equals the value `r`.
* A vector of length `D+1`. The first `D` elements of the
vector correspond to ``\\mathbf{a}`` while the last element is ``b``.
`PoincareMap` uses `ds`, higher order interpolation from DifferentialEquations.jl,
and root finding from Roots.jl, to create a high accuracy estimate of the section.
`PoincareMap` follows the [`DynamicalSystem`](@ref) interface with the following adjustments:
1. `dimension(pmap) == dimension(ds)`, even though the Poincaré
map is effectively 1 dimension less.
2. Like [`StroboscopicMap`](@ref) time is discrete and counts the iterations on the
surface of section. [`initial_time`](@ref) is always `0` and [`current_time`](@ref)
is current iteration number.
3. A new function [`current_crossing_time`](@ref) returns the real time corresponding
to the latest crossing of the hyperplane, which is what the [`current_state(ds)`](@ref)
corresponds to as well.
4. For the special case of `plane` being a `Tuple{Int, <:Real}`, a special `reinit!` method
is allowed with input state of length `D-1` instead of `D`, i.e., a reduced state already
on the hyperplane that is then converted into the `D` dimensional state.
## Example
```julia
using DynamicalSystemsBase
ds = Systems.rikitake(zeros(3); μ = 0.47, α = 1.0)
pmap = poincaremap(ds, (3, 0.0))
step!(pmap)
next_state_on_psos = current_state(pmap)
```
[^DatserisParlitz2022]:
Datseris & Parlitz 2022, _Nonlinear Dynamics: A Concise Introduction Interlaced with Code_,
[Springer Nature, Undergrad. Lect. Notes In Physics](https://doi.org/10.1007/978-3-030-91032-7)
"""
mutable struct PoincareMap{I<:ContinuousTimeDynamicalSystem, F, P, R, V} <: DiscreteTimeDynamicalSystem
ds::I
plane_distance::F
planecrossing::P
Tmax::Float64
rootkw::R
state_on_plane::V
tcross::Float64
t::Base.RefValue{Int}
# These two fields are for setting the state of the pmap from the plane
# (i.e., given a D-1 dimensional state, create the full D-dimensional state)
dummy::Vector{Float64}
diffidxs::Vector{Int}
end
Base.parent(pmap::PoincareMap) = pmap.ds
function PoincareMap(
ds::DS, plane;
Tmax = 1e3,
direction = -1, u0 = nothing,
rootkw = (xrtol = 1e-6, atol = 1e-8)
) where {DS<:ContinuousTimeDynamicalSystem}
reinit!(ds, u0)
D = dimension(ds)
check_hyperplane_match(plane, D)
planecrossing = PlaneCrossing(plane, direction > 0)
plane_distance = (t) -> planecrossing(ds(t))
v = recursivecopy(current_state(ds))
dummy = zeros(D)
diffidxs = _indices_on_poincare_plane(plane, D)
pmap = PoincareMap(
ds, plane_distance, planecrossing, Tmax, rootkw,
v, current_time(ds), Ref(0), dummy, diffidxs
)
step!(pmap)
pmap.t[] = 0 # first step is 0
return pmap
end
_indices_on_poincare_plane(plane::Tuple, D) = setdiff(1:D, [plane[1]])
_indices_on_poincare_plane(::Vector, D) = collect(1:D-1)
additional_details(pmap::PoincareMap) = [
"hyperplane" => pmap.planecrossing.plane,
"crossing time" => pmap.tcross,
]
###########################################################################################
# Extensions
###########################################################################################
for f in (:initial_state, :current_parameters, :initial_parameters, :referrenced_sciml_prob,
:dynamic_rule, :(SciMLBase.isinplace), :(StateSpaceSets.dimension))
@eval $(f)(pmap::PoincareMap, args...) = $(f)(pmap.ds, args...)
end
initial_time(pmap::PoincareMap) = 0
current_time(pmap::PoincareMap) = pmap.t[]
current_state(pmap::PoincareMap) = pmap.state_on_plane
"""
current_crossing_time(pmap::PoincareMap) → tcross
Return the time of the latest crossing of the Poincare section.
"""
current_crossing_time(pmap::PoincareMap) = pmap.tcross
function SciMLBase.step!(pmap::PoincareMap)
u, t = poincare_step!(pmap.ds, pmap.plane_distance, pmap.planecrossing, pmap.Tmax, pmap.rootkw)
if isnothing(u)
error("Exceeded `Tmax` without crossing the plane.")
else
pmap.state_on_plane = u # this is always a brand new vector
pmap.tcross = t
pmap.t[] = pmap.t[] + 1
return pmap
end
end
SciMLBase.step!(pmap::PoincareMap, n::Int, s = true) = (for _ ∈ 1:n; step!(pmap); end; pmap)
function SciMLBase.reinit!(pmap::PoincareMap, u0 = initial_state(pmap);
t0 = initial_time(pmap), p = current_parameters(pmap)
)
isnothing(u0) && return pmap
if length(u0) == dimension(pmap)
u = u0
elseif length(u0) == dimension(pmap) - 1
u = _recreate_state_from_poincare_plane(pmap, u0)
else
error("Dimension of state for `PoincareMap` reinit is inappropriate.")
end
reinit!(pmap.ds, u; t0, p)
step!(pmap) # always step once to reach the PSOS
pmap.t[] = 0 # first step is 0
pmap
end
function set_state!(pmap::PoincareMap, u)
set_state!(pmap.ds, u)
step!(pmap) # always step once to reach the PSOS
return
end
function _recreate_state_from_poincare_plane(pmap::PoincareMap, u0)
plane = pmap.planecrossing.plane
if plane isa Tuple
pmap.dummy[pmap.diffidxs] .= u0
pmap.dummy[plane[1]] = plane[2]
else
error("Don't know how to convert state on generic plane into full space.")
end
return pmap.dummy
end
###########################################################################################
# Poincare step
###########################################################################################
"""
poincare_step!(integ, plane_distance, planecrossing, Tmax, rootkw)
Low level function that actually performs the algorithm of finding the next crossing
of the Poincaré surface of section. Return the state and time at the section or `nothing` if
evolved for more than `Tmax` without any crossing.
"""
function poincare_step!(ds, plane_distance, planecrossing, Tmax, rootkw)
t0 = current_time(ds)
# Check if initial condition is already on the plane
side = planecrossing(current_state(ds))
if side == 0
dat = current_state(ds)
step!(ds)
return dat, t0
end
# Otherwise evolve until juuuuuust crossing the plane
tprev = t0
while side < 0
(current_time(ds) - t0) > Tmax && break
tprev = current_time(ds)
step!(ds)
side = planecrossing(current_state(ds))
end
while side ≥ 0
(current_time(ds) - t0) > Tmax && break
tprev = current_time(ds)
step!(ds)
side = planecrossing(current_state(ds))
end
# we evolved too long and no crossing, return nothing
(current_time(ds) - t0) > Tmax && return (nothing, nothing)
# Else, we're guaranteed to have time window before and after the plane
time_window = (tprev, current_time(ds))
tcross = Roots.find_zero(plane_distance, time_window, ROOTS_ALG; rootkw...)
ucross = ds(tcross)
return ucross, tcross
end
###########################################################################################
# Poincare surface of section
###########################################################################################
"""
poincaresos(ds::CoupledODEs, plane, T = 1000.0; kwargs...) → P::StateSpaceSet
Return the iterations of `ds` on the Poincaré surface of section with the `plane`,
by evolving `ds` up to a total of `T`.
Return a [`StateSpaceSet`](@ref) of the points that are on the surface of section.
This function initializes a [`PoincareMap`](@ref) and steps it until its
[`current_crossing_time`](@ref) exceeds `T`. You can also use [`trajectory`](@ref)
with [`PoincareMap`](@ref) to get a sequence of `N::Int` points instead.
The keywords `Ttr, save_idxs` act as in [`trajectory`](@ref).
See [`PoincareMap`](@ref) for `plane` and all other keywords.
"""
function poincaresos(ds::CoupledODEs, plane, T = 1000.0;
save_idxs = 1:dimension(ds), Ttr = 0, kwargs...
)
pmap = PoincareMap(ds, plane; kwargs...)
i = typeof(save_idxs) <: Int ? save_idxs : SVector{length(save_idxs), Int}(save_idxs...)
data = _initialize_output(current_state(pmap), i)
poincaresos!(data, pmap, i, T, Ttr)
end
function poincaresos!(data, pmap, i, T, Ttr)
tprev = current_crossing_time(pmap)
while current_crossing_time(pmap) - tprev < Ttr
step!(pmap)
end
push!(data, current_state(pmap))
tprev = current_crossing_time(pmap)
while current_crossing_time(pmap) - tprev < T
step!(pmap)
push!(data, current_state(pmap)[i])
end
return StateSpaceSet(data)
end
_initialize_output(::S, ::Int) where {S} = eltype(S)[]
_initialize_output(u::S, i::SVector{N, Int}) where {N, S} = typeof(u[i])[]
function _initialize_output(u, i)
error("The variable index when producing the PSOS must be Int or SVector{Int}")
end