Help with NonlinearSolve.jl

I am getting a cryptic error when trying to solve this system:

using StaticArrays
using NonlinearSolve

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

x0 = [0.1; 1.2]
x0s = MVector{size(x0)...}(x0)

prob = NonlinearProblem{true}(f!, x0s)
julia> solver = NonlinearSolve.solve(prob, NewtonRaphson(), tol=1e-9)
ERROR: UndefVarError: Tridiagonal not defined
Stacktrace:
 [1] (::NonlinearSolve.DefaultLinSolve)(x::MVector{2, Float64}, A::MMatrix{2, 2, Float64, 4}, b::MVector{2, Float64}, update_matrix::Bool; tol::Nothing, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ NonlinearSolve ~/.julia/packages/NonlinearSolve/KTnKV/src/utils.jl:151
 [2] DefaultLinSolve
   @ ~/.julia/packages/NonlinearSolve/KTnKV/src/utils.jl:133 [inlined]
 [3] perform_step(solver::NonlinearSolve.NewtonImmutableSolver{NonlinearFunction{true, typeof(f!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, NewtonRaphson{12, true, DataType, NonlinearSolve.DefaultLinSolve}, MVector{2, Float64}, MVector{2, Float64}, SciMLBase.NullParameters, typeof(NonlinearSolve.DEFAULT_NORM), Float64, NonlinearSolve.NewtonRaphsonCache{NonlinearSolve.JacobianWrapper{NonlinearFunction{true, typeof(f!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, SciMLBase.NullParameters}, NonlinearSolve.DefaultLinSolve, MMatrix{2, 2, Float64, 4}, MVector{2, Float64}, ForwardDiff.JacobianConfig{ForwardDiff.Tag{NonlinearSolve.JacobianWrapper{NonlinearFunction{true, typeof(f!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, SciMLBase.NullParameters}, Float64}, Float64, 2, Tuple{MVector{2, ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.JacobianWrapper{NonlinearFunction{true, typeof(f!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, SciMLBase.NullParameters}, Float64}, Float64, 2}}, MVector{2, ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.JacobianWrapper{NonlinearFunction{true, typeof(f!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, SciMLBase.NullParameters}, Float64}, Float64, 2}}}}}}, alg::NewtonRaphson{12, true, DataType, NonlinearSolve.DefaultLinSolve}, #unused#::Val{true})
   @ NonlinearSolve ~/.julia/packages/NonlinearSolve/KTnKV/src/raphson.jl:46
 [4] solve!(solver::NonlinearSolve.NewtonImmutableSolver{NonlinearFunction{true, typeof(f!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, NewtonRaphson{12, true, DataType, NonlinearSolve.DefaultLinSolve}, MVector{2, Float64}, MVector{2, Float64}, SciMLBase.NullParameters, typeof(NonlinearSolve.DEFAULT_NORM), Float64, NonlinearSolve.NewtonRaphsonCache{NonlinearSolve.JacobianWrapper{NonlinearFunction{true, typeof(f!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, SciMLBase.NullParameters}, NonlinearSolve.DefaultLinSolve, MMatrix{2, 2, Float64, 4}, MVector{2, Float64}, ForwardDiff.JacobianConfig{ForwardDiff.Tag{NonlinearSolve.JacobianWrapper{NonlinearFunction{true, typeof(f!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, SciMLBase.NullParameters}, Float64}, Float64, 2, Tuple{MVector{2, ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.JacobianWrapper{NonlinearFunction{true, typeof(f!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, SciMLBase.NullParameters}, Float64}, Float64, 2}}, MVector{2, ForwardDiff.Dual{ForwardDiff.Tag{NonlinearSolve.JacobianWrapper{NonlinearFunction{true, typeof(f!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, SciMLBase.NullParameters}, Float64}, Float64, 2}}}}}})
   @ NonlinearSolve ~/.julia/packages/NonlinearSolve/KTnKV/src/solve.jl:64
 [5] solve(::NonlinearProblem{MVector{2, Float64}, true, SciMLBase.NullParameters, NonlinearFunction{true, typeof(f!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::NewtonRaphson{12, true, DataType, NonlinearSolve.DefaultLinSolve}; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:tol,), Tuple{Float64}}})
   @ NonlinearSolve ~/.julia/packages/NonlinearSolve/KTnKV/src/solve.jl:5
 [6] top-level scope
   @ REPL[7]:1

I tried to follow the documentation, what am I missing? it seems that I need to load another package?

Open an issue. The error that you are getting is because the package’s nonlinear solvers weren’t setup to support MVectors at this time, but that can be fixed.

2 Likes

But then how do I work with functions that mutate their inputs, like the one above? If I use SArrays I get another error:

julia> solver = NonlinearSolve.solve(prob, NewtonRaphson(), tol=1e-9)
ERROR: setindex!(::SVector{2, Float64}, value, ::Int) is not defined.
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:33
 [2] setindex!(a::SVector{2, Float64}, value::Float64, i::Int64)
   @ StaticArrays ~/.julia/packages/StaticArrays/LJQEe/src/indexing.jl:3
 [3] f!(F::SVector{2, Float64}, x::SVector{2, Float64}, #unused#::SciMLBase.NullParameters)
   @ Main ./REPL[3]:2
 [4] (::NonlinearFunction{true, typeof(f!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing})(::SVector{2, Float64}, ::Vararg{Any})
   @ SciMLBase ~/.julia/packages/SciMLBase/DXiE6/src/scimlfunctions.jl:342
 [5] init(::NonlinearProblem{SVector{2, Float64}, true, SciMLBase.NullParameters, NonlinearFunction{true, typeof(f!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::NewtonRaphson{12, true, DataType, NonlinearSolve.DefaultLinSolve}; alias_u0::Bool, maxiters::Int64, tol::Float64, internalnorm::Function, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ NonlinearSolve ~/.julia/packages/NonlinearSolve/KTnKV/src/solve.jl:53
 [6] solve(::NonlinearProblem{SVector{2, Float64}, true, SciMLBase.NullParameters, NonlinearFunction{true, typeof(f!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::NewtonRaphson{12, true, DataType, NonlinearSolve.DefaultLinSolve}; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:tol,), Tuple{Float64}}})
   @ NonlinearSolve ~/.julia/packages/NonlinearSolve/KTnKV/src/solve.jl:4
 [7] top-level scope
   @ REPL[7]:1

Arrays work with mutation just fine.

I’ll fix the MAarrays issue rather soon though. The issue is just to get it in the email queue to be fixed.

1 Like

Yeah, I was trying to have it performant. I tested NLsolve vs NonlinearSolve and without static arrays, NLsolve is faster.

Than NonlinearSolve with SArray? That should be the fastest one since NLsolve with MArray will allocate buffers and such.

1 Like

NonlinearSolve with regular arrays vs NLsolve with regular arrays. NLsolve is faster in that case.

https://github.com/SciML/NonlinearSolve.jl/pull/50

using StaticArrays
using NonlinearSolve

function f(x, _)
    F1 = (x[1] + 3) * (x[2]^3 - 7) + 18
    F2 = sin(x[2] * exp(x[1]) - 1)
    SA[F1,F2]
end

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

x0 = [0.1; 1.2]
x0s = SVector{size(x0)...}(x0)
x0m = MVector{size(x0)...}(x0)
prob = NonlinearProblem{false}(f, x0s)

using NLsolve, BenchmarkTools
@btime sol = solve(prob,NewtonRaphson()); # 320.000 ns (2 allocations: 128 bytes)
@btime nlsolve(f!, x0m); # 1.460 μs (35 allocations: 1.36 KiB)

That is probably more of what you are looking for.

1 Like

Ohh, I didn’t try with that last line in f, I thought it would create an allocation, but it doesn’t.

What’s a good way to do the same as in f but with say 200 equations?

Arrays. Static arrays (both SVector and MVector) will not scale beyond about 10-20 equations.

1 Like

Right, but I would still need to feed NonlinearSolve.solve a non-mutating function, like f, instead of a mutating function like f!?
My question is how to best do that, even if it allocates. In other words, you use scalars F1, F2, how can I scale that to 200 equations?

EDIT:
Nevermind, I will just create F inside the function, populate it, then return it.