fixedpoints with automatic Jacobian (#311)

* fix problem with non given Jacobian

* bump version

CHANGELOG.md

@@ -1,5 +1,8 @@
 # main
 
+# v3.1
+- `fixedpoints` can now be called without a Jacobian (giving `nothing` as third argument, also the default value). In this case a Jacobian is estimated automatically via ForwardDiff.jl.
+
 # v3.0
 
 Major release part of the DynamicalSystems.jl v3.

Project.toml

@@ -1,7 +1,7 @@
 name = "ChaosTools"
 uuid = "608a59af-f2a3-5ad4-90b4-758bdf3122a7"
 repo = "https://github.com/JuliaDynamics/ChaosTools.jl.git"
-version = "3.0.2"
+version = "3.1.0"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"

src/stability/fixedpoints.jl

@@ -8,7 +8,7 @@ export fixedpoints, .., ×, IntervalBox, interval
 Return all fixed points `fp` of the given out-of-place `ds`
 (either `DeterministicIteratedMap` or `CoupledODEs`)
 that exist within the state space subset `box` for parameter configuration `p`.
-Fixed points are returned as a [`Dataset`](@ref).
+Fixed points are returned as a [`StateSpaceSet`](@ref).
 For convenience, a vector of the Jacobian eigenvalues of each fixed point, and whether
 the fixed points are stable or not, are also returned.
 
@@ -38,34 +38,33 @@ as a start of a continuation process. See also [`periodicorbits`](@ref).
 - `tol = 1e-15` is the root-finding tolerance.
 - `warn = true` throw a warning if no fixed points are found.
 """
-function fixedpoints(ds::DynamicalSystem, box, J;
+function fixedpoints(ds::DynamicalSystem, box, J = nothing;
         method = IntervalRootFinding.Krawczyk, tol = 1e-15, warn = true,
-        o = nothing, # the keyword `o` will be the period in a future version...
+        order = nothing, # the keyword `order` will be the period in a future version...
     )
     if isinplace(ds)
-        error("`fixedpoints` works only for out-of-place dynamical systems.")
+        error("`fixedpoints` currently works only for out-of-place dynamical systems.")
     end
     # Jacobian: copy code from `DynamicalSystemsBase`
     f = dynamic_rule(ds)
+    p = current_parameters(ds)
     if isnothing(J)
-        error("At the moment automatic differentiation doesn't work...")
-        Jf = (u, p, t) -> DynamicalSystemsBase.ForwardDiff.jacobian((x) -> f(x, p, t), u)
+        Jf(u, p, t) = DynamicalSystemsBase.ForwardDiff.jacobian(x -> f(x, p, 0.0), u)
     else
         Jf = J
     end
-    p = current_parameters(ds)
     # Find roots via IntervalRootFinding.jl
-    f = to_root_f(ds, p, o)
-    jac = to_root_J(Jf, ds, p, o)
-    r = IntervalRootFinding.roots(f, jac, box, method, tol)
+    fun = to_root_f(ds, p, order)
+    jac = to_root_J(Jf, ds, p, order)
+    r = IntervalRootFinding.roots(fun, jac, box, method, tol)
     D = dimension(ds)
-    fp::Dataset{D, Float64} = roots_to_dataset(r, D, warn)
+    fp = roots_to_dataset(r, D, warn)
     # Find eigenvalues and stability
     eigs = Vector{Vector{Complex{Float64}}}(undef, length(fp))
-    J = zeros(dimension(ds), dimension(ds)) # `eigvals` doesn't work with `SMatrix`
+    Jm = zeros(dimension(ds), dimension(ds)) # `eigvals` doesn't work with `SMatrix`
     for (i, u) in enumerate(fp)
-        J .= Jf(u, p, 0) # notice that we use the "pure" jacobian, no -u!
-        eigs[i] = LinearAlgebra.eigvals(Array(J))
+        Jm .= Jf(u, p, 0) # notice that we use the "pure" jacobian, no -u!
+        eigs[i] = LinearAlgebra.eigvals(Array(Jm))
     end
     stable = Bool[isstable(ds, e) for e in eigs]
     return fp, eigs, stable
@@ -98,7 +97,7 @@ end
 function roots_to_dataset(r, D, warn)
     if isempty(r) && warn
         @warn "No fixed points found!"
-        return Dataset{D, Float64}()
+        return StateSpaceSet{D, Float64}()
     end
     if any(root.status != :unique for root in r) && warn
         @warn "Non-unique fixed points found!"
@@ -107,7 +106,7 @@ function roots_to_dataset(r, D, warn)
     for (j, root) in enumerate(r)
         F[j, :] .= map(i -> (i.hi + i.lo)/2, root.interval)
     end
-    return Dataset(F; warn = false)
+    return StateSpaceSet(F; warn = false)
 end
 
 isstable(::CoupledODEs, e) = maximum(real(x) for x in e) < 0

test/stability/fixedpoints.jl

@@ -13,7 +13,6 @@ using ChaosTools, Test
     henon_fp = hmfp(current_parameters(ds)...)
 
     @testset "J=$(J)" for J in (nothing, henon_jacob)
-        isnothing(J) && continue
         fp, eigs, stable = fixedpoints(ds, box, J)
         @test size(fp) == (2,2)
         @test stable == [false, false]
@@ -75,7 +74,6 @@ end
     afps = lorenzfp(22,8/3)
 
     @testset "J=$(J)" for J in (nothing, lorenz_jacob)
-        isnothing(J) && continue
         fp, eigs, stable = fixedpoints(ds, box, J)
         @test length(fp) == 3
         for p in fp