Solving Non-linear Systems

It is often required to solve systems of non-linear equations in engineering. As for the the case of systems of linear equations, there are multiple options in Julia.

Available packages

A popular package is NLsolve.jl. This package provides three methods: Trust region, Newton’s method with optional line search, and Anderson acceleration. It can also solve fixed-point problems.

Alternatively, you could use JuMP - a very powerful mathematical programming modelling language, built on top of Julia, that can call several free and commercial solvers. In this case, set up a problem with non-linear constraints matching your system and a constant value objective function.

There are also several less popular packages available.

In the spirit of not wanting to learn several different interfaces to access options for different algorithms, we shall stick to the SciML package for non-linear systems, NonlinearSolve.jl.

This package uses the same interface as LinearSolve.jl. It has built-in in algorithms, but is also a unified front-end for other packages, including NLsolve.jl.

NonlinearSolve

The NonlinearSolve.jl package can be used to solve three types of problem:

  1. Interval root-finding (bracketing) problems, where a scalar root is found within a specified interval: find \(t ∈ [tₗ, tₕ]\) such that \(f(t) = 0\)
  2. Solving a system of non-linear equations, i.e. find \(\mathbf{u}\) such that \(f(\mathbf{u}) = \mathbf{0}\)
  3. Solving steady-state problems, i.e. solve \(u' = f(u, t) = f(u, \infty) = 0\)

To solve a problem, you need to provide a function that defines the system of equations. This can take one of two forms:

  • A non-mutating function (out-of-place), f(t, p), which returns u
  • A mutating function (in-place), f!(u, t, p), which mutates the parameter u

You then define either an IntervalNonlinearProblem, NonlinearProblem or SteadyStateProblem, depending on which of the three types of problems above you want to solve.

Finally, you call solve() to get the result.

As a simple demonstration, we shall calculate the boiling point of water, using the Antoine correlation with parameters from NIST. We calculate the boiling point by finding the root of a function that is the difference between the vapour pressure at a specified temperature and 1 atm. We do this in two ways, using the Newton-Raphson method from a starting guess and using the regula falsi method to find roots in a specified interval.

using NonlinearSolve

"""
    Calculate the difference between the vapour pressure and 1atm
"""
function f(u, p)
    A, B, C = p
    return 10^(A - B/(u + C)) - 1.01325
end

# Antoine parameters from NIST
# Stull, 1947
p = (4.6543, 1435.264, -64.848)

# Initial guess for NewtonRaphson
u0 = 373.15

# Specified span for Falsi
uspan = (370.0, 380.0);

prob = NonlinearProblem(f, u0, p)
sol = solve(prob, NewtonRaphson()) # u: 373.6009857758806

f(sol.u, p) # 4.121456953498637e-9

prob2 = IntervalNonlinearProblem(f, uspan, p)
sol2 = solve(prob2, Falsi()) # u: 373.60098565855054

f(sol2.u, p) # -1.5543122344752192e-15

We can also solve for the steady-state solution of a differential equation:

\[u'(t) = u(t) - u^2(t)\]

with \(u(0) = 1.5\). This ODE will have a steady-state at \(u = 1\).

using NonlinearSolve, DifferentialEquations

function ssode(u, p, t)
    return u - u^2
end

u0 = 1.5
probss = SteadyStateProblem(ssode, u0, nothing)
solss = solve(probss, DynamicSS(Vern6())) # u: 1.0000007379583993

While this is a simple demonstration of how to solve for an ODE’s steady-state, it would have been a lot more efficient to convert the problem to a non-linear equation, rather than using the approach of solving the dynamic ODE problem and running to steady-state.

Back to top