Add steady state docs and MINPACK

Author: ChrisRackauckas

docs/pages.jl

@@ -5,7 +5,9 @@ pages = ["index.md",
                        "tutorials/iterator_interface.md"],
     "Basics" => Any["basics/NonlinearProblem.md",
                     "basics/NonlinearFunctions.md",
+                    "basics/NonlinearSolution.md",
                     "basics/FAQ.md"],
     "Solvers" => Any["solvers/NonlinearSystemSolvers.md",
-                     "solvers/BracketingSolvers.md"],
+                     "solvers/BracketingSolvers.md",
+                     "solvers/SteadyStateSolvers.md"],
 ]

docs/src/basics/FAQ.md

@@ -35,4 +35,4 @@ MATLAB 2022a achieves 1.66s. Try this code yourself: we receive 0.06 seconds, or
 This example is still not optimized in the Julia code and we expect an improvement in a near
 future version.
 
-For more information on performance of SciML, see the [SciMLBenchmarks](https://docs.sciml.ai/SciMLBenchmarksOutput/stable/)
+For more information on performance of SciML, see the [SciMLBenchmarks](https://docs.sciml.ai/SciMLBenchmarksOutput/stable/).

docs/src/basics/NonlinearProblem.md

@@ -1,32 +1,46 @@
 # Nonlinear Problems
 
-## The Two Types of Nonlinear Problems
+## The Three Types of Nonlinear Problems
 
 NonlinearSolve.jl tackles two related types of nonlinear systems:
 
 1. Interval rootfinding problems. I.e., find the ``t in [t_0, t_f]`` such that `f(t) = 0`.
 2. Systems of nonlinear equations, i.e. find the `u` such that `f(u) = 0`.
+3. Steady state problems, i.e. find the `u` such that `u' = f(u,t)` has reached steady state,
+   i.e. `0 = f(u, ∞)`.
 
-The former is for solving scalar rootfinding problems, i.e. finding a single number, and
+The first is for solving scalar rootfinding problems, i.e. finding a single number, and
 requires that a bracketing interval is known. For a bracketing interval, one must have that
 the sign of `f(t_0)` is opposite the sign of `f(t_f)`, thus guaranteeing a root in the
 interval.
 
-The latter type of nonlinear system can be multidimensional and thus no ordering nor
-boundaries are assumed to be known. For a system of nonlinear equations, `f` can return
-an array and the solver seeks to find the value of `u` for which all outputs of `f` are
-simultaniously zero.
-
 !!! note
 
     Interval rootfinding problems allow for `f` to return an array, in which case the interval
     rootfinding problem is interpreted as finding the first `t` such that any of the components
     of the array hit zero.
 
+The second type of nonlinear system can be multidimensional and thus no ordering nor
+boundaries are assumed to be known. For a system of nonlinear equations, `f` can return
+an array and the solver seeks to find the value of `u` for which all outputs of `f` are
+simultaniously zero.
+
+The last type if equivalent to a nonlinear system but with the extra interpretion of
+having a potentially preferred unique root. That is, when there are multiple `u` such
+that `f(u) = 0`, the `NonlinearProblem` does not have a preferred solution, while for the
+`SteadyStateProblem` the preferred solution is the `u(∞)` that would arise from solving the
+ODE `u' = f(u,t)`.
+
+!!! warn
+
+    Most solvers for `SteadyStateProblem` do not guarentee the preferred solution and
+    instead will solve for some `u` in the set of solutions. The documentation of the
+    nonlinear solvers will note if they return the preferred solution.
 
 ## Problem Construction Details
 
 ```@docs
 SciMLBase.IntervalNonlinearProblem
 SciMLBase.NonlinearProblem
+SciMLBase.SteadyStateProblem
 ```

docs/src/basics/NonlinearSolution.md

@@ -0,0 +1,5 @@
+# Nonlinear Solutions
+
+```@docs
+SciMLBase.NonlinearSolution
+```

docs/src/solvers/NonlinearSystemSolvers.md

@@ -17,7 +17,7 @@ is an implementation which is specialized for small equations. It is non-allocat
 static arrays and thus really well-optimized for small systems, thus usually outperforming
 the other methods when such types are used for `u0`. `NLSolveJL`'s
 `:trust_region` method can be a good choice for high stability, along with
-`CMINPACK`.s
+`DynamicSS`.
 
 For a system which is very non-stiff (i.e., the condition number of the Jacobian
 is small, or the eigenvalues of the Jacobian are within a few orders of magnitude),
@@ -51,9 +51,48 @@ can be used directly to reduce dependencies and improve load times.
   very efficient and non-allocating. Thus this implmentation is well-suited for small
   systems of equations.
 
+### SteadyStateDiffEq.jl
+
+Note that these solvers do not come by default, and thus one needs to install
+the package before using these solvers:
+
+```julia
+]add SteadyStateDiffEq
+using SteadyStateDiffEq
+```
+
+SteadyStateDiffEq.jl uses ODE solvers to iteratively approach the steady state. It is a
+very stable method for solving nonlinear systems, though in many cases can be more
+computationally expensive than direct methods.
+
+- `DynamicSS` : Uses an ODE solver to find the steady state. Automatically
+  terminates when close to the steady state.
+  `DynamicSS(alg;abstol=1e-8,reltol=1e-6,tspan=Inf)` requires that an
+  ODE algorithm is given as the first argument.  The absolute and
+  relative tolerances specify the termination conditions on the
+  derivative's closeness to zero.  This internally uses the
+  `TerminateSteadyState` callback from the Callback Library.  The
+  simulated time for which given ODE is solved can be limited by
+  `tspan`.  If `tspan` is a number, it is equivalent to passing
+  `(zero(tspan), tspan)`.
+
+Example usage:
+
+```julia
+using NonlinearSolve, SteadyStateDiffEq, OrdinaryDiffEq
+sol = solve(prob,DynamicSS(Tsit5()))
+
+using Sundials
+sol = solve(prob,DynamicSS(CVODE_BDF()),dt=1.0)
+```
+
+!!! note
+
+    If you use `CVODE_BDF` you may need to give a starting `dt` via `dt=....`.*
+
 ### SciMLNLSolve.jl
 
-This is a wrapper package for importing solvers from other packages into this interface.
+This is a wrapper package for importing solvers from NLsolve.jl into this interface.
 Note that these solvers do not come by default, and thus one needs to install
 the package before using these solvers:
 
@@ -62,7 +101,6 @@ the package before using these solvers:
 using SciMLNLSolve
 ```
 
-- `CMINPACK()`: A wrapper for using the classic MINPACK method through [MINPACK.jl](https://github.com/sglyon/MINPACK.jl)
 - `NLSolveJL()`: A wrapper for [NLsolve.jl](https://github.com/JuliaNLSolvers/NLsolve.jl)
 
 ```julia
@@ -134,3 +172,34 @@ The choices for the linear solver are:
 - `:KLU`: A sparse factorization method. Requires that the user specify a
   Jacobian. The Jacobian must be set as a sparse matrix in the `ODEProblem`
   type.
+
+### NonlinearSolveMINPACK
+
+This is a wrapper package for importing solvers from MINPACK.jl into this interface.
+Note that these solvers do not come by default, and thus one needs to install
+the package before using these solvers:
+
+```julia
+]add NonlinearSolveMINPACK
+using NonlinearSolveMINPACK
+```
+
+- `CMINPACK()`: A wrapper for using the classic MINPACK method through [MINPACK.jl](https://github.com/sglyon/MINPACK.jl)
+
+```julia
+CMINPACK(;show_trace::Bool=false, tracing::Bool=false, method::Symbol=:hybr,
+          io::IO=stdout)
+```
+
+The keyword argument `method` can take on different value depending on which method of `fsolve` you are calling. The standard choices of `method` are:
+
+- `:hybr`: Modified version of Powell's algorithm. Uses MINPACK routine [`hybrd1`](https://github.com/devernay/cminpack/blob/d1f5f5a273862ca1bbcf58394e4ac060d9e22c76/hybrd1.c)
+- `:lm`: Levenberg-Marquardt. Uses MINPACK routine [`lmdif1`](https://github.com/devernay/cminpack/blob/d1f5f5a273862ca1bbcf58394e4ac060d9e22c76/lmdif1.c)
+- `:lmdif`: Advanced Levenberg-Marquardt (more options available with `;kwargs...`). See MINPACK routine [`lmdif`](https://github.com/devernay/cminpack/blob/d1f5f5a273862ca1bbcf58394e4ac060d9e22c76/lmdif.c) for more information
+- `:hybrd`: Advacned modified version of Powell's algorithm (more options available with `;kwargs...`). See MINPACK routine [`hybrd`](https://github.com/devernay/cminpack/blob/d1f5f5a273862ca1bbcf58394e4ac060d9e22c76/hybrd.c) for more information
+
+If a Jacobian is supplied as part of the [`NonlinearFunction`](@ref nonlinearfunctions),
+then the following methods are allowed:
+
+- `:hybr`: Advacned modified version of Powell's algorithm with user supplied Jacobian. Additional arguments are available via `;kwargs...`. See MINPACK routine [`hybrj`](https://github.com/devernay/cminpack/blob/d1f5f5a273862ca1bbcf58394e4ac060d9e22c76/hybrj.c) for more information
+- `:lm`: Advanced Levenberg-Marquardt with user supplied Jacobian. Additional arguments are available via `;kwargs...`. See MINPACK routine [`lmder`](https://github.com/devernay/cminpack/blob/d1f5f5a273862ca1bbcf58394e4ac060d9e22c76/lmder.c) for more information

docs/src/solvers/SteadyStateSolvers.md

@@ -0,0 +1,62 @@
+# Steady State Solvers
+
+`solve(prob::SteadyStateProblem,alg;kwargs)`
+
+Solves for the steady states in the problem defined by `prob` using the algorithm
+`alg`. If no algorithm is given, a default algorithm will be chosen.
+
+## Recommended Methods
+
+Conversion to a NonlinearProblem is generally the fastest method. However, this will not
+guarantee the preferred root, and thus if the preferred root is required, then it's
+recommended that one uses `DynamicSS`. For `DynamicSS`, in many cases an adaptive stiff
+solver, like a Rosenbrock or BDF method (`Rodas5` or `QNDF`), is a good way to allow for
+very large time steps as the steady state approaches.
+
+!!! note
+
+    The SteadyStateDiffEq.jl methods on a `SteadyStateProblem` respect the time definition
+    in the nonlinear definition, i.e. `u' = f(u,t)` uses the correct values for `t` as the
+    solution evolves. A conversion of a `SteadyStateProblem` to a `NonlinearProblem`
+    replaces this with the nonlinear system `u' = f(u,∞)`, and thus the direct
+    `SteadyStateProblem` approach can give different answers (i.e. the correct unique
+    fixed point) on ODEs with non-autonomous dynamics.
+
+## Full List of Methods
+
+### Conversion to NonlinearProblem
+
+Any `SteadyStateProblem` can be trivially converted to a `NonlinearProblem` via
+`NonlinearProblem(prob::SteadyStateProblem)`. Using this appraoch, any of the solvers from
+the [Nonlinear System Solvers page](@ref nonlinearsystemsolvers) can be used.
+
+### SteadyStateDiffEq.jl
+
+SteadyStateDiffEq.jl uses ODE solvers to iteratively approach the steady state. It is a
+very stable method for solving nonlinear systems, though in many cases can be more
+computationally expensive than direct methods.
+
+- `DynamicSS` : Uses an ODE solver to find the steady state. Automatically
+  terminates when close to the steady state.
+  `DynamicSS(alg;abstol=1e-8,reltol=1e-6,tspan=Inf)` requires that an
+  ODE algorithm is given as the first argument.  The absolute and
+  relative tolerances specify the termination conditions on the
+  derivative's closeness to zero.  This internally uses the
+  `TerminateSteadyState` callback from the Callback Library.  The
+  simulated time for which given ODE is solved can be limited by
+  `tspan`.  If `tspan` is a number, it is equivalent to passing
+  `(zero(tspan), tspan)`.
+
+Example usage:
+
+```julia
+using NonlinearSolve, SteadyStateDiffEq, OrdinaryDiffEq
+sol = solve(prob,DynamicSS(Tsit5()))
+
+using Sundials
+sol = solve(prob,DynamicSS(CVODE_BDF()),dt=1.0)
+```
+
+!!! note
+
+    If you use `CVODE_BDF` you may need to give a starting `dt` via `dt=....`.*