Support for fixed points of n-th order in fixedpoints (#326)

* support fixed points of n-th iterate of a map

* test

* simplified jacobian

* test fix

src/stability/fixedpoints.jl

@@ -37,14 +37,29 @@ as a start of a continuation process. See also [`periodicorbits`](@ref).
   see the docs of IntervalRootFinding.jl for all possibilities.
 - `tol = 1e-15` is the root-finding tolerance.
 - `warn = true` throw a warning if no fixed points are found.
+- `order = nothing` search for fixed points of the n-th iterate of 
+  [`DeterministicIteratedMap`](@ref). Must be a positive integer or `nothing`.
+  Select `nothing` or 1 to search for the fixed points of the original map.
+
+## Performance notes
+
+Setting `order` to a value greater than 5 can be very slow. Consider using 
+more suitable algorithms for periodic orbit detection, such as 
+[`periodicorbits`](@ref).
+
 """
 function fixedpoints(ds::DynamicalSystem, box, J = nothing;
         method = IntervalRootFinding.Krawczyk, tol = 1e-15, warn = true,
-        order = nothing, # the keyword `order` will be the period in a future version...
+        order = nothing,
     )
     if isinplace(ds)
         error("`fixedpoints` currently works only for out-of-place dynamical systems.")
     end
+
+    if !(isnothing(order) || (isa(order, Int) && order > 0))
+        error("`order` must be a positive integer or `nothing`.")
+    end
+
     # Jacobian: copy code from `DynamicalSystemsBase`
     f = dynamic_rule(ds)
     p = current_parameters(ds)
@@ -83,15 +98,28 @@ end
 function to_root_f(ds::DeterministicIteratedMap, p, o::Int)
     f = dynamic_rule(ds)
     u -> begin
-        v = copy(u) # copy is free for `SVector`
+        v = u
         for _ in 1:o
             v = f(v, p, 0)
         end
         return v - u
     end
 end
-# TODO: Estimate periodic orbits of period `o` for discrete systems.
-# What's left is to create the Jacobian for higher iterates.
+
+# this can be derived from chain rule : J(f(g(x))) = Jf(g(x))*Jg(x)
+function to_root_J(Jf, ds::DeterministicIteratedMap, p, o::Int)
+    c = Diagonal(ones(typeof(current_state(ds))))
+    d = dimension(ds)
+    u -> begin
+        J = Jf(u, p, 0.0)
+        x = ds.f(u, p, 0.0)
+        for i in 2:o
+            J = Jf(x, p, 0.0) * J
+            x = ds.f(x, p, 0.0)
+        end
+        return J .- c
+    end
+end
 
 
 function roots_to_dataset(r, D, warn)

test/stability/fixedpoints.jl

@@ -38,6 +38,23 @@ end
 #         end
 #     end
 # end
+@testset "Henon map (n-th order)" begin
+    henon_rule(x, p, n) = SVector{2}(1.0 - p[1]*x[1]^2 + x[2], p[2]*x[1])
+    henon_jacob(x, p, n) = SMatrix{2,2}(-2*p[1]*x[1], p[2], 1.0, 0.0)
+    ds = DeterministicIteratedMap(henon_rule, zeros(2), [1.4, 0.3])
+    x = interval(-3.0, 3.0)
+    y = interval(-10.0, 10.0)
+    box = x × y
+    o = 4
+    fp, eigs, stable = fixedpoints(ds, box, henon_jacob; order=o)
+    tol = 1e-8
+    for x0 in fp
+        set_state!(ds, x0)
+        step!(ds, o)
+        xn = current_state(ds)
+        @test isapprox(x0, xn; atol = tol)
+    end
+end
 
 @testset "Lorenz system" begin
     function lorenz_rule(u, p, t)