|
| 1 | +# Getting Started with Nonlinear Rootfinding in Julia |
| 2 | + |
| 3 | +NonlinearSolve.jl is a system for solving rootfinding problems, i.e. finding |
| 4 | +the value $$u$$ such that $$f(u) = 0$$. In this tutorial we will go through |
| 5 | +the basics of NonlinearSolve.jl, demonstrating the core ideas and leading you |
| 6 | +to understanding the deeper parts of the documentation. |
| 7 | + |
| 8 | +## The Three Types of Nonlinear Systems |
| 9 | + |
| 10 | +There are three types of nonlinear systems: |
| 11 | + |
| 12 | +1. The "standard nonlinear system", i.e. the `NonlinearProblem`. This is a |
| 13 | + system of equations with an initial condition where you want to satisfy |
| 14 | + all equations simultaniously. |
| 15 | +2. The "interval rootfinding problem", i.e. the `IntervalNonlinearProblem`. |
| 16 | + This is the case where you're given an interval `[a,b]` and need to find |
| 17 | + where `f(u) = 0` for `u` inside the bounds. |
| 18 | +3. The "steady state problem", i.e. find the `u` such that `u' = f(u) = 0`. |
| 19 | + While related to (1), it's not entirely the same because there's a uniquely |
| 20 | + defined privledged root. |
| 21 | +4. The nonlinear least squares problem, which is an overconstrained nonlinear |
| 22 | + system (i.e. more equations than states) which might not be satisfiable, i.e. |
| 23 | + there may be no `u` such that `f(u) = 0`, and thus we find the `u` which |
| 24 | + minimizes `||f(u)||` in the least squares sense. |
| 25 | + |
| 26 | +For now let's focus on the first two. The other two are covered in later tutorials, |
| 27 | +but from the first two we can show the general flow of the NonlinearSolve.jl package. |
| 28 | + |
| 29 | +## Problem Type 1: Solving Nonlinear Systems of Equations |
| 30 | + |
| 31 | +A nonlinear system $$f(u) = 0$$ is specified by defining a function `f(u,p)`, |
| 32 | +where `p` are the parameters of the system. For example, the following solves |
| 33 | +the vector equation $$f(u) = u^2 - p$$ for a vector of equations: |
| 34 | + |
| 35 | +```@example |
| 36 | +using NonlinearSolve |
| 37 | +
|
| 38 | +f(u, p) = u .* u .- p |
| 39 | +u0 = [1.0, 1.0] |
| 40 | +p = 2.0 |
| 41 | +prob = NonlinearProblem(f, u0, p) |
| 42 | +sol = solve(prob) |
| 43 | +``` |
| 44 | + |
| 45 | +where `u0` is the initial condition for the rootfinder. Native NonlinearSolve.jl |
| 46 | +solvers use the given type of `u0` to determine the type used within the solver |
| 47 | +and the return. Note that the parameters `p` can be any type, but most are an |
| 48 | +AbstractArray for automatic differentiation. |
| 49 | + |
| 50 | +### Investigating the Solution |
| 51 | + |
| 52 | +To investigate the solution, one can look at the elements of the `NonlinearSolution`. |
| 53 | +The most important value is `sol.u`: this is the `u` that satisfies `f(u) = 0`. For example: |
| 54 | + |
| 55 | +```@example |
| 56 | +u = sol.u |
| 57 | +``` |
| 58 | + |
| 59 | +```@example |
| 60 | +f(u, p) |
| 61 | +``` |
| 62 | + |
| 63 | +This final value, the difference of the solution against zero, can also be found with `sol.resid`: |
| 64 | + |
| 65 | +```@example |
| 66 | +sol.resid |
| 67 | +``` |
| 68 | + |
| 69 | +To know if the solution converged, or why the solution had not converged we can check the return |
| 70 | +code (`retcode`): |
| 71 | + |
| 72 | +```@example |
| 73 | +sol.retcode |
| 74 | +``` |
| 75 | + |
| 76 | +There are multiple return codes which can mean the solve was successful, and thus we can use the |
| 77 | +general command `SciMLBase.successful_retcode` to check whether the solution process exited as |
| 78 | +intended: |
| 79 | + |
| 80 | +```@exmaple |
| 81 | +SciMLBase.successful_retcode(sol) |
| 82 | +``` |
| 83 | + |
| 84 | +If we're curious about what it took to solve this equation, then you're in luck because all of the |
| 85 | +details can be found in `sol.stats`: |
| 86 | + |
| 87 | +```@example |
| 88 | +sol.stats |
| 89 | +``` |
| 90 | + |
| 91 | +For more information on `NonlinearSolution`s, see the [`NonlinearSolution` manual page](@ref solution). |
| 92 | + |
| 93 | +### Interacting with the Solver Options |
| 94 | + |
| 95 | +While `sol = solve(prob)` worked for our case here, in many situations you may need to interact more |
| 96 | +deeply with how the solving process is done. First things first, you can specify the solver using the |
| 97 | +positional arguments. For example, let's set the solver to `TrustRegion()`: |
| 98 | + |
| 99 | +```@example |
| 100 | +solve(prob, TrustRegion()) |
| 101 | +``` |
| 102 | + |
| 103 | +For a complete list of solver choices, see [the nonlinear system solvers page](@ref nonlinearsystemsolvers). |
| 104 | + |
| 105 | +Next we can modify the tolerances. Here let's set some really low tolerances to force a tight solution: |
| 106 | + |
| 107 | +```@example |
| 108 | +solve(prob, TrustRegion(), reltol=1e-12, abstol=1e-12) |
| 109 | +``` |
| 110 | + |
| 111 | +There are many more options for doing this configuring. Specifically for handling termination conditions, |
| 112 | +see the [Termination Conditions](@ref termination_condition) page for more details. And for more details on |
| 113 | +all of the available keyword arguments, see the [solver options](@ref solver_options) page. |
| 114 | + |
| 115 | +## Problem Type 2: Solving Interval Rootfinding Problems with Bracketing Methods |
| 116 | + |
| 117 | +For scalar rootfinding problems, bracketing methods exist in NonlinearSolve. The difference with bracketing |
| 118 | +methods is that with bracketing methods, instead of giving a `u0` initial condition, you pass a `uspan (a,b)` |
| 119 | +bracket in which the zero is expected to live. For example: |
| 120 | + |
| 121 | +```@example |
| 122 | +using NonlinearSolve |
| 123 | +f(u, p) = u * u - 2.0 |
| 124 | +uspan = (1.0, 2.0) # brackets |
| 125 | +prob_int = IntervalNonlinearProblem(f, uspan) |
| 126 | +sol = solve(prob_int) |
| 127 | +``` |
| 128 | + |
| 129 | +All of the same option handling from before works just as before, now just with different solver choices |
| 130 | +(see the [bracketing solvers](@ref bracketing) page for more details). For example, let's set the solver |
| 131 | +to `ITP()` and set a high absolute tolerance: |
| 132 | + |
| 133 | +```@example |
| 134 | +sol = solve(prob_int, ITP(), abstol = 0.01) |
| 135 | +``` |
| 136 | + |
| 137 | +## Going Beyond the Basics: How to use the Documentation |
| 138 | + |
| 139 | +Congrats, you now know how to use the basics of NonlinearSolve.jl! However, there is so much more to |
| 140 | +see. Next check out: |
| 141 | + |
| 142 | +- [Some code optimization tricks to know about with NonlinearSolve.jl](@ref code_optimization) |
| 143 | +- [An iterator interface which lets you step through the solving process step by step](@ref iterator) |
| 144 | +- [How to handle large systems of equations efficiently](@ref large_systems) |
| 145 | +- [Ways to use NonlinearSolve.jl that is faster to startup and can statically compile](@ref fast_startup) |
| 146 | +- [More detailed termination conditions](@ref termination_conditions_tutorial) |
| 147 | + |
| 148 | +And also check out the rest of the manual. |
0 commit comments