Rank-deficient jacobian error in Newton-Raphson method for unstable limit cycle refinement

Specific Domains Numerics
, ,
1

I am trying to refine the position of an unstable limit cycle within a chaotic attractor using the Newton-Raphson method, but I am running into a confusing error related to a rank-deficient matrix.

To avoid degeneracy of the Jacobian matrix, at each iteration I project the initial point onto a secant plane. This secant plane is chosen as the direction corresponding to the largest absolute component of the phase velocity vector. This technique should remove local non-uniqueness of the solution and prevent the Jacobian from becoming singular. (The details of the secant plane choice are clear from the code.)

In practice, for a pretty good initial guess of a periodic orbit in the Lorenz system

Begin state = [-1.9598067770796535, 2.0076976037190195, 26.58579174986873]
End state   = [-1.9578169398768586, 2.002913423511019, 26.580466855969398]
Δ = 0.007429826875196052, T = 39.44

the method reports a rank-deficient matrix error:

Algorithm: NewtonRaphson(
    descent = NewtonDescent(),
    autodiff = AutoForwardDiff(),
    vjp_autodiff = AutoFiniteDiff(
        fdtype = Val{:forward}(),
        fdjtype = Val{:forward}(),
        fdhtype = Val{:hcentral}(),
        dir = true
    ),
    jvp_autodiff = AutoForwardDiff(),
    concrete_jac = Val{false}()
)

----            -------------           -----------             -------
Iter            f(u) inf-norm           Step 2-norm             cond(J)
----            -------------           -----------             -------
0               8.30327189e-312         0.00000000e+00          Inf
┌ Warning: Potential Rank Deficient Matrix Detected. Attempting to solve using Pivoted QR Factorization.
└
1               8.30327189e-312         0.00000000e+00          Inf
Final           8.30327189e-312
----------------------

From a mathematical point of view, I did not expect any problems to arise.

I would appreciate any help!
Below is the MWE I am using:

using NonlinearSolve
using DynamicalSystems
using ForwardDiff
using LinearAlgebra


function newton_method(ds::CoupledODEs, u0::U where {U<:AbstractArray{<:Real}}, T::Float64, Δt::Float64)
    # Dimensionality of the system
    D = dimension(ds)

    # First of all we'll solve the problem for systems of dimension 3
    @assert(D == 3)

    # v - the value of the phase velocity vector at point u0
    v = ds.integ.f.f(u0, ds.p0, 0.0)

    # k - index of the component that determines the position of the secant
    k = argmax(abs.(v))

    # i and j - indices of the remaining components that will change
    i, j = filter((x) -> (return x != k), [i for i in 1:D])[1:2]

    # Inital state
    x0 = [u0[[i, j]]; T]

    # Detefinition of the residual
    f = (err, x, p) -> begin
      # Creating a vector of initial conditions
      u = SVector{D}(
        if k == 1
          (u0[k], x[1], x[2])
        elseif k == 2
          (x[1], u0[k], x[2])
        elseif k == 3
          (x[1], x[2], u0[k])
        end
      )

      # Estimated time of the period
      T_local = x[end]

      # Ranges in which we will integrate the system
      tspan = (0.0, T_local)

      # Integration of the system for the initial condition u for time T_local
      sol = solve(remake(ds.integ.sol.prob; u0 = u, tspan = tspan);
                  DynamicalSystemsBase.DEFAULT_DIFFEQ..., ds.diffeq...)

      # The difference between the final and initial state
      err = sol.u[end] .- sol.u[1]
    end

    prob = NonlinearProblem(f, x0, nothing)

    sol = solve(prob, NewtonRaphson();
                show_trace = Val(true), trace_level = TraceAll(10))

    u0_new = zeros(3)
    u0_new[i] = sol.u[1]
    u0_new[j] = sol.u[2]
    u0_new[k] = u0[k]

    T_new  = sol.u[end]

    return (states = trajectory(ds, T_new - Δt, u0_new; Δt = Δt)[1], period = T_new)
end

function main()
  system = Systems.lorenz()

  # Initial guess
  ig = [-1.9598067770796535, 2.0076976037190195, 26.58579174986873]
  # Period guess
  T  = 39.44
  Δt = 0.01

  res = trajectory(system, T - Δt, ig; Δt = Δt)[1]
  Δ   = LinearAlgebra.norm(res[end] .- res[1])

  println("Begin state = $ig")
  println("End state   = $(res[end])")
  println("Δ = $Δ, T = $T")

  res_newton = newton_method(system, ig, T, Δt)

  Δ_newton = LinearAlgebra.norm(res_newton.states[end] .- res_newton.states[1])
  println("Begin state after newton method = $(res_newton.states[1])")
  println("Δ_newton = $Δ_newton, T = $(res_newton.period)")
end

main()
2

Why are you choosing this method? This is exactly the situation where you don’t use NewtonRaphson and instead use something like a TrustRegion or line search to stabilize it. What does the default method give?

1 Like
3

From my experience simple standard shooting does not work well on unstable orbits, You should use multiple shooting. This is avaible in BifurcationKit.jl. See for example the tutorials

2 Likes
4

There is this implementation of an algorithm developed for finding periodic orbits in chaotic attractors:

it is based on nonlinear solve.

1 Like