# Rank-deficient jacobian error in Newton-Raphson method for unstable limit cycle refinement
**URL:** https://discourse.julialang.org/t/rank-deficient-jacobian-error-in-newton-raphson-method-for-unstable-limit-cycle-refinement/129975
**Category:** Numerics
**Tags:** newton-raphson, nonlinearsolve, dynamical-systems
**Created:** [June 17, 2025, 6:14pm UTC](https://discourse.julialang.org/t/rank-deficient-jacobian-error-in-newton-raphson-method-for-unstable-limit-cycle-refinement/129975 "2025-06-17T18:14:41Z")
**Posts on this page:** 4
**Page:** 1
### Author:  [@R.Kanakov](https://discourse.julialang.org/u/R.Kanakov)
#### Post date: [June 17, 2025, 6:14pm UTC](https://discourse.julialang.org/t/rank-deficient-jacobian-error-in-newton-raphson-method-for-unstable-limit-cycle-refinement/129975/1 "2025-06-17T18:14:41Z")
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
```plaintext
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:
```plaintext
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:
```julia
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()
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [June 17, 2025, 7:24pm UTC](https://discourse.julialang.org/t/rank-deficient-jacobian-error-in-newton-raphson-method-for-unstable-limit-cycle-refinement/129975/2 "2025-06-17T19:24:02Z")
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?
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [June 17, 2025, 9:56pm UTC](https://discourse.julialang.org/t/rank-deficient-jacobian-error-in-newton-raphson-method-for-unstable-limit-cycle-refinement/129975/3 "2025-06-17T21:56:36Z")
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](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/tutorials/tutorials/)
---
### Author:  [@Datseris](https://discourse.julialang.org/u/Datseris)
#### Post date: [June 17, 2025, 10:20pm UTC](https://discourse.julialang.org/t/rank-deficient-jacobian-error-in-newton-raphson-method-for-unstable-limit-cycle-refinement/129975/4 "2025-06-17T22:20:57Z")
There is this implementation of an algorithm developed for finding periodic orbits in chaotic attractors:
> **[Tutorial · PeriodicOrbits.jl](https://juliadynamics.github.io/PeriodicOrbits.jl/dev/tutorial/)**
>
> Documentation for PeriodicOrbits.jl.
it is based on nonlinear solve.