Modeling Nonlinear Systems

ModelingToolkit.jl is not only useful for generating initial value problems (ODEProblem). The package can also build nonlinear systems. This is, for example, useful for finding the steady state of an ODE. This steady state is reached when the nonlinear system of differential equations equals zero.

Note

The high level @mtkmodel macro used in the getting started tutorial is not yet compatible with NonlinearSystem. We thus have to use a lower level interface to define nonlinear systems. For an introduction to this interface, read the programmatically generating Systems tutorial.

using ModelingToolkit, NonlinearSolve

# Define a nonlinear system
@variables x y z
@parameters σ ρ β
eqs = [0 ~ σ * (y - x)
       0 ~ x * (ρ - z) - y
       0 ~ x * y - β * z]
@mtkcompile ns = System(eqs)

guesses = [x => 1.0, y => 0.0, z => 0.0]
ps = [σ => 10.0, ρ => 26.0, β => 8 / 3]

prob = NonlinearProblem(ns, vcat(guesses, ps))
sol = solve(prob, NewtonRaphson())

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

We found the x, y and z for which the right hand sides of eqs are all equal to zero.

Just like with ODEProblems we can generate the NonlinearProblem with its analytical Jacobian function:

prob = NonlinearProblem(ns, vcat(guesses, ps), jac = true)
sol = solve(prob, NewtonRaphson())

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