Solving Nonlinear Systems with NLsolve

In this example we will go one step deeper and showcase the direct function generation capabilities in ModelingToolkit.jl to build nonlinear systems. Let's say we wanted to solve for the steady state of the previous ODE. This is the nonlinear system defined by where the derivatives are zero. We use (unknown) variables for our nonlinear system.

using ModelingToolkit

@variables x y z
@parameters σ ρ β

# Define a nonlinear system
eqs = [0 ~ σ*(y-x),
       0 ~ x*(ρ-z)-y,
       0 ~ x*y - β*z]
ns = NonlinearSystem(eqs, [x,y,z], [σ,ρ,β])
nlsys_func = generate_function(ns)[2] # second is the inplace version

which generates:

(var"##MTIIPVar#405", u, p)->begin
        @inbounds begin
                @inbounds begin
                        let (x, y, z, σ, ρ, β) = (u[1], u[2], u[3], p[1], p[2], p[3])
                            var"##MTIIPVar#405"[1] = (*)(σ, (-)(y, x))
                            var"##MTIIPVar#405"[2] = (-)((*)(x, (-)(ρ, z)), y)
                            var"##MTIIPVar#405"[3] = (-)((*)(x, y), (*)(β, z))
                        end
                    end
            end
        nothing
    end

We can use this to build a nonlinear function for use with NLsolve.jl:

f = eval(nlsys_func)
du = zeros(3); u = ones(3)
params = (10.0,26.0,2.33)
f(du,u,params)
du

#=
3-element Array{Float64,1}:
  0.0
 24.0
 -1.33
 =#

We can similarly ask to generate the in-place Jacobian function:

j_func = generate_jacobian(ns)[2] # second is in-place
j! = eval(j_func)

which gives:

:((var"##MTIIPVar#582", u, p)->begin
          #= C:\Users\accou\.julia\dev\ModelingToolkit\src\utils.jl:70 =#
          #= C:\Users\accou\.julia\dev\ModelingToolkit\src\utils.jl:71 =#
          #= C:\Users\accou\.julia\dev\ModelingToolkit\src\utils.jl:71 =# @inbounds begin
                  #= C:\Users\accou\.julia\dev\ModelingToolkit\src\utils.jl:72 =#
                  #= C:\Users\accou\.julia\dev\ModelingToolkit\src\utils.jl:53 =# @inbounds begin
                          #= C:\Users\accou\.julia\dev\ModelingToolkit\src\utils.jl:53 =#
                          let (x, y, z, σ, ρ, β) = (u[1], u[2], u[3], p[1], p[2], p[3])
                              var"##MTIIPVar#582"[1] = (*)(σ, -1)
                              var"##MTIIPVar#582"[2] = (-)(ρ, z)
                              var"##MTIIPVar#582"[3] = y
                              var"##MTIIPVar#582"[4] = σ
                              var"##MTIIPVar#582"[5] = -1
                              var"##MTIIPVar#582"[6] = x
                              var"##MTIIPVar#582"[7] = 0
                              var"##MTIIPVar#582"[8] = (*)(x, -1)
                              var"##MTIIPVar#582"[9] = (*)(-1, β)
                          end
                      end
              end
          #= C:\Users\accou\.julia\dev\ModelingToolkit\src\utils.jl:74 =#
          nothing
      end)

Now, we can call nlsolve by enclosing our parameters into the functions:

using NLsolve
nlsolve((out, x) -> f(out, x, params), (out, x) -> j!(out, x, params), ones(3))

If one would like the generated function to be a Julia function instead of an expression, and allow this function to be used from within the same world-age, one simply needs to pass Val{false} to tell it to generate the function, i.e.:

nlsys_func = generate_function(ns, [x,y,z], [σ,ρ,β], expression=Val{false})[2]

which uses GeneralizedGenerated.jl to build the same world-age function on the fly without eval.