Getting Started

Suppose the system of nonlinear equations of interest can be described as follows:

using NonlinearSystems

# Residual function
function f!(F, x)
    F[1] = (x[1] + 3) * (x[2]^3 - 7) + 18
    F[2] = sin(x[2] * exp(x[1]) - 1)
    return F
end

# Jacobian function (optional)
function j!(J, x)
    J[1,1] = x[2]^3 - 7
    J[1,2] = 3 * x[2]^2 * (x[1] + 3)
    u = exp(x[1]) * cos(x[2] * exp(x[1]) - 1)
    J[2,1] = x[2] * u
    J[2,2] = u
    return J
end

# Initial value
x0 = [0.1, 1.2]

To solve the above equations as a root-finding problem, we specify the Hybrid algorithm by passing either Hybrid or Hybrid{RootFinding} as the first argument:

# Evaluate Jacobians via finite differencing methods from FiniteDiff.jl
solve(Hybrid, f!, x0)
# Use user-specified Jacobian function and separate out the initialization step
s = init(Hybrid, f!, j!, x0)
solve!(s)

The last line from above calls a non-allocating method solve! that mutates the pre-allocated problem s in-place. On Julia REPL, the essential information is summarized as follows:

julia> s2×2 NonlinearSystem{RootFinding, Vector{Float64}, Matrix{Float64}, HybridSolver{Float64, DenseLUSolver{LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}}, Vector{Float64}}, Nothing, Nothing}:
  Problem type:                 Root finding
  Algorithm:                    Hybrid
  Candidate (x):                [-6.06188e-10, 1.0]
  ‖f(x)‖₂:                      6.971672272635292e-9
  ‖Ddx‖₂:                       1.1197356805172518e-5
  Solver exit state:            ftol_reached
  Iterations:                   7
  Residual calls (f):           7
  Jacobian calls (df/dx):       1

The solution can be retrieved by accessing the corresponding field:

julia> s.x2-element Vector{Float64}:
 -6.061879300097132e-10
  0.9999999988463306

In general, a nonlinear system of equations may have multiple solutions or no solution. If only solutions that fall in certain regions are of interest, we may impose lower and upper bounds that force the solver to only seek solution candidates within the bounded areas:

julia> solve(Hybrid, f!, ones(2), lower=fill(0.5,2), upper=fill(2.0,2))2×2 NonlinearSystem{RootFinding, Vector{Float64}, Matrix{Float64}, HybridSolver{Float64, DenseLUSolver{LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}}, Vector{Float64}}, Vector{Float64}, Vector{Float64}}:
  Problem type:                 Root finding
  Algorithm:                    Hybrid
  Candidate (x):                [1.10117, 1.37701]
  ‖f(x)‖₂:                      1.66896278574161e-9
  ‖Ddx‖₂:                       1.0356454220559748e-6
  Solver exit state:            ftol_reached
  Iterations:                   8
  Residual calls (f):           11
  Jacobian calls (df/dx):       3

Notice that we have found a different solution in the above example.

In practice, we often do not expect the existence of solutions that solve the equations exactly as identities. Instead, we may solve a least-squares problem that minimizes the Euclidean norm of residuals. To do so, simply specify the algorithm type as Hybrid{LeastSquares}:

julia> s = solve(Hybrid{LeastSquares}, f!, x0)2×2 NonlinearSystem{LeastSquares, Vector{Float64}, Matrix{Float64}, HybridSolver{Float64, DenseCholeskySolver{Float64, Int8, Matrix{Float64}, Vector{Float64}, Nothing}, Vector{Float64}}, Nothing, Nothing}:
  Problem type:                 Least squares
  Algorithm:                    Hybrid
  Candidate (x):                [-6.06188e-10, 1.0]
  ‖f(x)‖₂:                      6.971668834975478e-9
  ‖∇f(x)‖₂:                     NaN
  ‖Ddx‖₂:                       1.119735680119611e-5
  Solver exit state:            ftol_reached
  Iterations:                   7
  Residual calls (f):           7
  Jacobian calls (df/dx):       1

Notice that the gradient norm is NaN. For this specific problem, convergence is attained before the gradient is ever evaluated.

To inspect the solver iteration, summary information can be printed for each evaluation:

julia> s = solve(Hybrid{LeastSquares}, f!, x0, showtrace=1);  iter    1  =>  ‖f(x)‖₂ =  1.687505e+00  ‖∇f(x)‖₂ =           NaN  ‖Ddx‖₂ =           NaN  δ =  1.612854e+01  ρ =           NaN
  iter    2  =>  ‖f(x)‖₂ =  3.755760e-01  ‖∇f(x)‖₂ =           NaN  ‖Ddx‖₂ =  2.350344e+00  δ =  2.350344e+00  ρ =  7.774372e-01
  iter    3  =>  ‖f(x)‖₂ =  5.511315e-02  ‖∇f(x)‖₂ =           NaN  ‖Ddx‖₂ =  3.312116e-01  δ =  4.700688e+00  ρ =  8.532570e-01
  iter    4  =>  ‖f(x)‖₂ =  2.152558e-03  ‖∇f(x)‖₂ =           NaN  ‖Ddx‖₂ =  3.953811e-02  δ =  4.700688e+00  ρ =  9.609429e-01
  iter    5  =>  ‖f(x)‖₂ =  5.613648e-04  ‖∇f(x)‖₂ =           NaN  ‖Ddx‖₂ =  4.914400e-03  δ =  7.907623e-02  ρ =  7.392103e-01
  iter    6  =>  ‖f(x)‖₂ =  6.682932e-06  ‖∇f(x)‖₂ =           NaN  ‖Ddx‖₂ =  9.781638e-04  δ =  7.907623e-02  ρ =  9.880952e-01

Notice that all relevant information is collected in a single object:

NonlinearSystem(::Type{P}, fdf::OnceDifferentiable, x, fx, dx, solver; kwargs...)

Instead of calling solve or solve!, which simply iterates NonlinearSystem in a loop, we may manually iterate the solver steps as follows:

julia> s = init(Hybrid{LeastSquares}, f!, x0);
julia> s.solverHybridSolver{Float64, DenseCholeskySolver{Float64, Int8, Matrix{Float64}, Vector{Float64}, Nothing}, Vector{Float64}}:
  iter    1  =>  ‖f(x)‖₂ =  1.687505e+00  ‖∇f(x)‖₂ =           NaN  ‖Ddx‖₂ =           NaN  δ =  1.612854e+01  ρ =           NaN
julia> iterate(s)(2×2 NonlinearSystem{LeastSquares}(Hybrid, 0.3755760088573572, inprogress), (2, NonlinearSystems.normal, NonlinearSystems.inprogress))
julia> s.solverHybridSolver{Float64, DenseCholeskySolver{Float64, Int8, Matrix{Float64}, Vector{Float64}, Nothing}, Vector{Float64}}:
  iter    2  =>  ‖f(x)‖₂ =  3.755760e-01  ‖∇f(x)‖₂ =           NaN  ‖Ddx‖₂ =  2.350344e+00  δ =  2.350344e+00  ρ =  7.774372e-01