# 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 `SArray`s 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.

``````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.