SciML / NonlinearSolve.jl

High-performance and differentiation-enabled nonlinear solvers (Newton methods), bracketed rootfinding (bisection, Falsi), with sparsity and Newton-Krylov support.
https://docs.sciml.ai/NonlinearSolve/stable/
MIT License
235 stars 41 forks source link

Storing trace #58

Closed KronosTheLate closed 11 months ago

KronosTheLate commented 2 years ago

I want to store the trace during solving. I would love to use the native NonlinearSolve methods for that, but it appears to me that no such thing is implemented. So I want to use the NLSolveJL interface. But from reading the documentation, I am still not able to see - how do i use it?

So far I have tried the following

#The functional version:
julia> solver = solve(probN, NewtonRaphson(), tol = 1e-9)
u: 8-element Vector{Float64}:
  27.52327682597094
  27.523962383296873
   3.2064613688401113
  13.258386042046885
   0.8000852068277563
   0.4578205135984268
  27.929656903088112
 124.54137097709595

#Failed attempts:
julia> solver = solve(probN, NLSolveJL(method=:newton,), tol = 1e-9)
ERROR: UndefVarError: NLSolveJL not defined
Stacktrace:
 [1] top-level scope
   @ REPL[9]:1

julia> solver = NLSolveJL(probN, tol = 1e-9, method=:newton)        
ERROR: UndefVarError: NLSolveJL not defined
Stacktrace:
 [1] top-level scope
   @ REPL[11]:1

With NLSolveJL not being defined, I was wondering if it is simply not exported. But this seems not to be the case:

julia> NonlinearSolve.NLSolveJL
ERROR: UndefVarError: NLSolveJL not defined
Stacktrace:
 [1] getproperty(x::Module, f::Symbol)
   @ Base .\Base.jl:35
 [2] top-level scope
   @ REPL[12]:1

I think that a section with working examples in the docs of how to access non-native solvers would go a long ways.

ChrisRackauckas commented 2 years ago

Did you do using NLsolve?

KronosTheLate commented 2 years ago

In a fresh session:

julia> using NonlinearSolve

julia>

julia> v = 120  #kg
120

julia> k = 2.5  #m
2.5

julia> w = 4.0  #km/m
4.0

julia> α = 2e-7 #kg⁻¹
2.0e-7

julia>

julia> function F(x⃗, parameters)
           L₀, L, p, x, θ, φ, a, H = tuple(x⃗...)
           return [
           a*(cosh(x/a)-1) - p,
           2a * sinh(x/a) - L ,
           2x + 2k*cos(θ) - parameters.d,
           p+k*sin(θ) - parameters.n,
           sinh(x/a) - tan(φ),
           (1+v/(w*L₀)) * tan(φ) - tan(θ),
           L₀ * (1 + α*H) - L,
           w*L₀ / 2sin(φ) - H,
           ]
       end
F (generic function with 1 method)

julia> params = (d=30, n=5)
(d = 30, n = 5)

julia> x₀ =    [29, 27, 1, params.d/2, deg2rad(45),  deg2rad(22.5), 40, 500]
8-element Vector{Float64}:
  29.0
  27.0
   1.0
  15.0
   0.7853981633974483
   0.39269908169872414
  40.0
 500.0

julia> probN = NonlinearProblem{false}(F, x₀, params)
NonlinearProblem with uType Vector{Float64}. In-place: false
u0: 8-element Vector{Float64}:
  29.0
  27.0
   1.0
  15.0
   0.7853981633974483
   0.39269908169872414
  40.0
 500.0

julia> solve(probN, NewtonRaphson(), tol = 1e-9)
u: 8-element Vector{Float64}:
  27.52327682597094
  27.523962383296873
   3.2064613688401113
  13.258386042046885
   0.8000852068277563
   0.4578205135984268
  27.929656903088112
 124.54137097709595

julia>

julia> using NLsolve

julia> solve(probN, NLSolveJL(method = :newton,), tol = 1e-9)
ERROR: UndefVarError: NLSolveJL not defined
Stacktrace:
 [1] top-level scope
   @ REPL[12]:1

julia> NLSolveJL(probN, method = :newton, tol = 1e-9)
ERROR: UndefVarError: NLSolveJL not defined
Stacktrace:
 [1] top-level scope
   @ REPL[13]:1
ChrisRackauckas commented 2 years ago

https://github.com/SciML/SciMLNLSolve.jl

It looks like that is documented? http://nonlinearsolve.sciml.ai/dev/solvers/NonlinearSystemSolvers/#SciMLNLSolve.jl

KronosTheLate commented 2 years ago

It was - I was not carefull enough in reading the documentation. The full working example is now

julia> using NonlinearSolve

julia>

julia> v = 120  #kg
120

julia> k = 2.5  #m
2.5

julia> w = 4.0  #km/m
4.0

julia> α = 2e-7 #kg⁻¹
2.0e-7

julia>

julia> function F(x⃗, parameters)
           L₀, L, p, x, θ, φ, a, H = tuple(x⃗...)
           return [
           a*(cosh(x/a)-1) - p,
           2a * sinh(x/a) - L ,
           2x + 2k*cos(θ) - parameters.d,
           p+k*sin(θ) - parameters.n,
           sinh(x/a) - tan(φ),
           (1+v/(w*L₀)) * tan(φ) - tan(θ),
           L₀ * (1 + α*H) - L,
           w*L₀ / 2sin(φ) - H,
           ]
       end
F (generic function with 1 method)

julia> params = (d=30, n=5)
(d = 30, n = 5)

julia> x₀ =    [29, 27, 1, params.d/2, deg2rad(45),  deg2rad(22.5), 40, 500]
8-element Vector{Float64}:
  29.0
  27.0
   1.0
  15.0
   0.7853981633974483
   0.39269908169872414
  40.0
 500.0

julia> probN = NonlinearProblem{false}(F, x₀, params)
NonlinearProblem with uType Vector{Float64}. In-place: false
u0: 8-element Vector{Float64}:
  29.0
  27.0
   1.0
  15.0
   0.7853981633974483
   0.39269908169872414
  40.0
 500.0

julia> using SciMLNLSolve

julia> test = solve(probN, NLSolveJL(method = :newton, show_trace=true), tol = 1e-9)
Iter     f(x) inf-norm    Step 2-norm
------   --------------   --------------
     0     3.484387e+02              NaN
     1     3.846384e+00     1.430607e+05
     2     2.695128e-02     1.122826e+01
     3     2.280210e-04     5.025951e-02
     4     1.677076e-08     3.750412e-06
u: 8-element Vector{Float64}:
  27.52327682125824
  27.523962378584063
   3.2064613689525356
  13.258386041929604
   0.8000852067627725
   0.4578205135111699
  27.929656760847312
 124.54137097770474

I think that for other careless readers like me, a simple example of how to use other solvers would be instructive.

ChrisRackauckas commented 2 years ago

Yup, we can add a tutorial

avik-pal commented 11 months ago
using NonlinearSolve

v = 120  #kg

k = 2.5  #m

w = 4.0  #km/m

α = 2e-7 #kg⁻¹

function F(x⃗, parameters)
    L₀, L, p, x, θ, φ, a, H = tuple(x⃗...)
    return [
        a * (cosh(x / a) - 1) - p,
        2a * sinh(x / a) - L,
        2x + 2k * cos(θ) - parameters.d,
        p + k * sin(θ) - parameters.n,
        sinh(x / a) - tan(φ),
        (1 + v / (w * L₀)) * tan(φ) - tan(θ),
        L₀ * (1 + α * H) - L,
        w * L₀ / 2sin(φ) - H,
    ]
end

params = (d = 30, n = 5)
x₀ = [29, 27, 1, params.d / 2, deg2rad(45), deg2rad(22.5), 40, 500]

probN = NonlinearProblem{false}(F, x₀, params)

sol = solve(probN; abstol = 1e-9, show_trace = Val(true), trace_level = TraceAll(100),
    store_trace = Val(true))

sol.trace will hold the trace of the solve process and it should display something like

Algorithm: GeneralKlement(linsolve = LUFactorization())

----     -------------        -----------          -------             
Iter     f(u) inf-norm        Step 2-norm          cond(J)             
----     -------------        -----------          -------             
0        3.48438696e+02       0.00000000e+00       1.00000000e+00      
Final    4.58673147e+00      
----------------------      

Algorithm: GeneralBroyden()

----     -------------        -----------          -------             
Iter     f(u) inf-norm        Step 2-norm          cond(J)             
----     -------------        -----------          -------             
0        3.48438696e+02       0.00000000e+00       1.00000000e+00      
100      2.21817021e+02       3.67955669e+01       8.24919384e+03      
Final    5.63586054e+00      
----------------------      

Algorithm: NewtonRaphson(ad = AutoForwardDiff())

----     -------------        -----------          -------             
Iter     f(u) inf-norm        Step 2-norm          cond(J)             
----     -------------        -----------          -------             
0        3.48438696e+02       0.00000000e+00       Inf                 
Final    3.55271368e-15      
----------------------      
u: 8-element Vector{Float64}:
  27.52327682597094
  27.523962383296873
   3.2064613688401113
  13.258386042046885
   0.8000852068277563
   0.4578205135984268
  27.929656903088112
 124.54137097709595