Help with NonlinearSolve.jl

General Usage
,
1

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?

2

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
3

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
4

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
5

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

6

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

1 Like
7

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

8

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
9

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

10

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

11

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

1 Like
12

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.