Find the root of an equation (i.e. solve f(u)=0)

A nonlinear system $f(u) = 0$ is specified by defining a function f(u,p), where p are the parameters of the system. Many problems can be written in such a way that solving a nonlinear rootfinding problem gives the solution. For example:

• Do you want to know $u$ such that $4^u + 6^u = 7^u$? Then solve $f(u) = 4^u + 6^u - 7^u = 0$ for u!
• If you have an ODE $u' = f(u)$, what is the point where the solution will be completely still, i.e. u' = 0?

All of these problems are solved by using a numerical rootfinder. Let's solve our first rootfind problem!

Required Dependencies

The following parts of the SciML Ecosystem will be used in this tutorial:

Problem Setup

For example, the following solves the vector equation:

\[\begin{aligned} 0 &= σ*(y-x)\\ 0 &= x*(ρ-z)-y\\ 0 &= x*y - β*z\\ \end{aligned}\]

With the parameter values $\sigma = 10.0$, $\rho = 26.0$, $\beta = 8/3$.

# Import the packages
using ModelingToolkit, NonlinearSolve

# Define the nonlinear system
@variables x=1.0 y=0.0 z=0.0
@parameters σ=10.0 ρ=26.0 β=8 / 3

eqs = [0 ~ σ * (y - x),
    0 ~ x * (ρ - z) - y,
    0 ~ x * y - β * z]
@mtkbuild ns = NonlinearSystem(eqs, [x, y, z], [σ, ρ, β])

# Convert the symbolic system into a numerical system
prob = NonlinearProblem(ns, [])

# Solve the numerical problem
sol = solve(prob, NewtonRaphson())

# Analyze the solution
@show sol[[x, y, z]], sol.resid

([-2.129924444096732e-29, -5.537803554651503e-28, -2.398137151871876e-28], [-5.32481111024183e-27, 6.395032404991669e-28])

Step-by-Step Solution

Step 1: Import the Packages

To do this tutorial, we will need a few components:

• ModelingToolkit.jl, our modeling environment
• NonlinearSolve.jl, the nonlinear system solvers

To start, let's add these packages as demonstrated in the installation tutorial:

using Pkg
Pkg.add(["ModelingToolkit", "NonlinearSolve"])

Now we're ready. Let's load in these packages:

# Import the packages
using ModelingToolkit, NonlinearSolve

Step 2: Define the Nonlinear System

Now let's define our nonlinear system. We use the ModelingToolkit.@variabes statement to declare our 3 state variables:

# Define the nonlinear system
@variables x=1.0 y=0.0 z=0.0

\[ \begin{equation} \left[ \begin{array}{c} x \\ y \\ z \\ \end{array} \right] \end{equation} \]

Notice that we are using the form state = initial condition. This is a nice shorthand for coupling an initial condition to our states. We now must similarly define our parameters, which we can associate default values via the form parameter = default value. This looks like:

@parameters σ=10.0 ρ=26.0 β=8 / 3

\[ \begin{equation} \left[ \begin{array}{c} \sigma \\ \rho \\ \beta \\ \end{array} \right] \end{equation} \]

Now we create an array of equations to define our nonlinear system that must be satisfied. This looks as follows:

eqs = [0 ~ σ * (y - x),
    0 ~ x * (ρ - z) - y,
    0 ~ x * y - β * z]

\[ \begin{align} 0 &= \left( - x + y \right) \sigma \\ 0 &= - y + x \left( - z + \rho \right) \\ 0 &= x y - z \beta \end{align} \]

Finally, we bring these pieces together, the equation along with its states and parameters, define our NonlinearSystem:

@mtkbuild ns = NonlinearSystem(eqs, [x, y, z], [σ, ρ, β])

\[ \begin{align} 0 &= \left( - x + y \right) \sigma \\ 0 &= x y - z \beta \end{align} \]

Step 3: Convert the Symbolic Problem to a Numerical Problem

Now that we have created our system, let's turn it into a numerical problem to approximate. This is done with the NonlinearProblem constructor, that transforms it from a symbolic ModelingToolkit representation to a numerical NonlinearSolve representation. We need to tell it the numerical details for whether to override any of the default values for the initial conditions and parameters.

In this case, we will use the default values for all our variables, so we will pass a blank override []. This looks like:

# Convert the symbolic system into a numerical system
prob = NonlinearProblem(ns, [])

NonlinearProblem with uType Vector{Float64}. In-place: true
u0: 2-element Vector{Float64}:
 1.0
 0.0

If we did want to change the initial condition of x to 2.0 and the parameter σ to 4.0, we would do [x => 2.0, σ => 4.0]. This looks like:

prob2 = NonlinearProblem(ns, [x => 2.0, σ => 4.0])

NonlinearProblem with uType Vector{Float64}. In-place: true
u0: 2-element Vector{Float64}:
 2.0
 0.0

Step 4: Solve the Numerical Problem

Now we solve the nonlinear system. For this, we choose a solver from the NonlinearSolve.jl's solver options. We will choose NewtonRaphson as follows:

# Solve the numerical problem
sol = solve(prob, NewtonRaphson())

retcode: Success
u: 2-element Vector{Float64}:
 -2.129924444096732e-29
 -2.398137151871876e-28

Step 5: Analyze the Solution

Now let's check out the solution. First of all, what kind of thing is the sol? We can see that by asking for its type:

typeof(sol)

SciMLBase.NonlinearSolution{Float64, 1, Vector{Float64}, Vector{Float64}, SciMLBase.NonlinearProblem{Vector{Float64}, true, ModelingToolkit.MTKParameters{Vector{Float64}, Vector{Float64}, Tuple{}, Tuple{}, Tuple{}, Tuple{}}, SciMLBase.NonlinearFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.GeneratedFunctionWrapper{(2, 2, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x7fc35376, 0xffe1a926, 0xf77ab257, 0x8fa63362, 0xe64c72b9), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf22f4034, 0x209d6fb0, 0x8c7b27bc, 0x03ebc810, 0x64336046), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ModelingToolkit.NonlinearSystem}, Nothing, ModelingToolkit.NonlinearSystem, Nothing, SciMLBase.OverrideInitData{SciMLBase.NonlinearProblem{Nothing, true, ModelingToolkit.MTKParameters{Vector{Float64}, StaticArraysCore.SizedVector{0, Float64, Vector{Float64}}, Tuple{}, Tuple{}, Tuple{}, Tuple{}}, SciMLBase.NonlinearFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.GeneratedFunctionWrapper{(2, 2, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xeda91b60, 0x8c51d865, 0xbf94bb2b, 0x632e456b, 0xbad8322a), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x1a82c6c9, 0x1e286c09, 0x152497a0, 0x4fde61ed, 0x654d8d84), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ModelingToolkit.NonlinearSystem}, Nothing, ModelingToolkit.NonlinearSystem, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, Nothing, Nothing, Nothing, ModelingToolkit.InitializationMetadata{ModelingToolkit.ReconstructInitializeprob{ModelingToolkit.var"#_getter#802"{Tuple{ComposedFunction{ModelingToolkit.PConstructorApplicator{typeof(identity)}, ComposedFunction{Base.Fix1{typeof(reduce), typeof(vcat)}, SymbolicIndexingInterface.MultipleGetters{Nothing, Tuple{SymbolicIndexingInterface.MultipleGetters{Nothing, Vector{SymbolicIndexingInterface.GetParameterIndex{ModelingToolkit.ParameterIndex{SciMLStructures.Tunable, Int64}}}}, SymbolicIndexingInterface.MultipleGetters{Nothing, Vector{SymbolicIndexingInterface.GetParameterIndex{ModelingToolkit.ParameterIndex{SciMLStructures.Initials, Int64}}}}}}}}, Returns{StaticArraysCore.SizedVector{0, Float64, Vector{Float64}}}, Returns{Tuple{}}, Returns{Tuple{}}, Returns{Tuple{}}}}, ComposedFunction{typeof(identity), SymbolicIndexingInterface.MultipleGetters{SymbolicIndexingInterface.ContinuousTimeseries, Vector{Any}}}}, Nothing, ModelingToolkit.SetInitialUnknowns{SymbolicIndexingInterface.MultipleSetters{Vector{SymbolicIndexingInterface.ParameterHookWrapper{SymbolicIndexingInterface.SetParameterIndex{ModelingToolkit.ParameterIndex{SciMLStructures.Initials, Int64}}, SymbolicUtils.BasicSymbolic{Real}}}}}}, Val{true}}}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithm{Missing, Missing, NonlinearSolveBase.NewtonDescent{Nothing}, ADTypes.AutoForwardDiff{nothing, Nothing}, ADTypes.AutoFiniteDiff{Val{:forward}, Val{:forward}, Val{:hcentral}, Nothing, Nothing, Bool}, ADTypes.AutoForwardDiff{nothing, Nothing}, Val{false}}, Nothing, Nothing, SciMLBase.NLStats, NonlinearSolveBase.NonlinearSolveTrace{Val{false}, Val{false}, Nothing, NonlinearSolveBase.NonlinearSolveTracing{Val{:minimal}}, SciMLBase.NonlinearProblem{Vector{Float64}, true, ModelingToolkit.MTKParameters{Vector{Float64}, Vector{Float64}, Tuple{}, Tuple{}, Tuple{}, Tuple{}}, SciMLBase.NonlinearFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.GeneratedFunctionWrapper{(2, 2, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x7fc35376, 0xffe1a926, 0xf77ab257, 0x8fa63362, 0xe64c72b9), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf22f4034, 0x209d6fb0, 0x8c7b27bc, 0x03ebc810, 0x64336046), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ModelingToolkit.NonlinearSystem}, Nothing, ModelingToolkit.NonlinearSystem, Nothing, SciMLBase.OverrideInitData{SciMLBase.NonlinearProblem{Nothing, true, ModelingToolkit.MTKParameters{Vector{Float64}, StaticArraysCore.SizedVector{0, Float64, Vector{Float64}}, Tuple{}, Tuple{}, Tuple{}, Tuple{}}, SciMLBase.NonlinearFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.GeneratedFunctionWrapper{(2, 2, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xeda91b60, 0x8c51d865, 0xbf94bb2b, 0x632e456b, 0xbad8322a), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x1a82c6c9, 0x1e286c09, 0x152497a0, 0x4fde61ed, 0x654d8d84), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{ModelingToolkit.NonlinearSystem}, Nothing, ModelingToolkit.NonlinearSystem, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, Nothing, Nothing, Nothing, ModelingToolkit.InitializationMetadata{ModelingToolkit.ReconstructInitializeprob{ModelingToolkit.var"#_getter#802"{Tuple{ComposedFunction{ModelingToolkit.PConstructorApplicator{typeof(identity)}, ComposedFunction{Base.Fix1{typeof(reduce), typeof(vcat)}, SymbolicIndexingInterface.MultipleGetters{Nothing, Tuple{SymbolicIndexingInterface.MultipleGetters{Nothing, Vector{SymbolicIndexingInterface.GetParameterIndex{ModelingToolkit.ParameterIndex{SciMLStructures.Tunable, Int64}}}}, SymbolicIndexingInterface.MultipleGetters{Nothing, Vector{SymbolicIndexingInterface.GetParameterIndex{ModelingToolkit.ParameterIndex{SciMLStructures.Initials, Int64}}}}}}}}, Returns{StaticArraysCore.SizedVector{0, Float64, Vector{Float64}}}, Returns{Tuple{}}, Returns{Tuple{}}, Returns{Tuple{}}}}, ComposedFunction{typeof(identity), SymbolicIndexingInterface.MultipleGetters{SymbolicIndexingInterface.ContinuousTimeseries, Vector{Any}}}}, Nothing, ModelingToolkit.SetInitialUnknowns{SymbolicIndexingInterface.MultipleSetters{Vector{SymbolicIndexingInterface.ParameterHookWrapper{SymbolicIndexingInterface.SetParameterIndex{ModelingToolkit.ParameterIndex{SciMLStructures.Initials, Int64}}, SymbolicUtils.BasicSymbolic{Real}}}}}}, Val{true}}}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}}}

From this, we can see that it is an NonlinearSolution. We can see the documentation for how to use the NonlinearSolution by checking the NonlinearSolve.jl solution type page. For example, the solution is stored as .u. What is the solution to our nonlinear system, and what is the final residual value? We can check it as follows:

# Analyze the solution
@show sol[[x, y, z]], sol.resid

([-2.129924444096732e-29, -5.537803554651503e-28, -2.398137151871876e-28], [-5.32481111024183e-27, 6.395032404991669e-28])