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()