NonlinearEquations

module ModiaMath.NonlinearEquations

This module contains functions to solve nonlinear algebraic equations:

• solveOneNonlinearEquation: Function that computes the solution of one non-linear algebraic equation y=f(u) in one unknown u in a reliable and efficient way using Brents algorithm.

• KINSOL: Module containing functions to solve a system of nonlinear systems of equations with Sundials KINSOL. The module is designed so that the same system is solved several times with KINSOL as needed by Modia simulations (auxiliary memory is only allocated once and not for every call). KINSOL is used in ModiaMath to solve nonlinear algebraic equations during initialization and at events of a simulation.

Main developer

Martin Otter, DLR - Institute of System Dynamics and Control

y = solveOneNonlinearEquation(f, u_min, u_max; u_nominal=1.0, tolerance=100.0*eps())

This function determines the solution of one non-linear algebraic equation y=f(u) in one unknown u in a reliable and efficient way. It is one of the best numerical algorithms for this purpose. As input, the nonlinear function f(u) has to be given, as well as an interval u_min, u_max that contains the solution, that is f(u_min) and f(u_max) must have a different sign. The function computes iteratively smaller intervals in which a sign change is present and terminates when the following test is fulfilled (u_nominal is the nominal value of u, so the order of magnitude of u around the solution and tolerance is the relative tolerance; for example, tolerance=1e-10 means that the solution is correct up to 10 digits):

absTol = 0.1*u_nominal*tolerance
if abs(length_of_u_interval) <= max(tolerance*abs(u), absTol) || f(u) == 0.0
    # root u found (interval is small enough)

The interval reduction is performed using inverse quadratic interpolation (interpolating with a quadratic polynomial through the last three points and computing the zero). If this fails, bisection is used, which always reduces the interval by a factor of two. The inverse quadratic interpolation method has superlinear convergence. This is roughly the same convergence rate as a globally convergent Newton method, but without the need to compute derivatives of the non-linear function. The solver function is a direct mapping of the Algol 60 procedure zero to Julia, from:

• Brent R.P. (1973): Algorithms for Minimization without derivatives. Prentice Hall, 1973, pp. 58-59.
  Download: http://wwwmaths.anu.edu.au/~brent/pd/rpb011i.pdf
  Errata and new print: http://wwwmaths.anu.edu.au/~brent/pub/pub011.html

Examples

import ModiaMath: solveOneNonlinearEquation

fun(u; w=1.0) = 3*u - sin(w*u) - 1
u_zero = solveOneNonlinearEquation(u->fun(u; w=3.0), 0.0, 5.0)
println("residue = ", fun(u_zero; w=3.0))

Remarks

• The interface was made such that it is identical to function Modelica.Math.Nonlinear.solveOneNonlinearEquation in order that automatic translation of Modelica to Modia is simplified. However, the termination condition was changed: The original Brent algorithm uses an absolute tolerance for termination (abs(lengthofinterval) <= 2eps()abs(u) + tolerance || f(u) == 0.0), whereas this was changed here to take a relative tolerance into account and use a similar definition of tolerances as used by Modia and Modelica, because easier to understand by a user (with relative tolerance and nominal value, instead of relative and absolute tolerances).

• Newer algorithms for the problem are presented in Alefeld, Potra, Shi (1995): Algorithm 748: enclosing zeros of continuous functions. Here, an inverse cubic interpolation is used instead of an inverse quadratic interpolation as in case of the Brent algorithm. The numerical experiments in this article with 15 test problems show that Brents algorithm needs about 5% more function evaluations for all the 15 test problems as the newly presented algorithm 4.2 (in some cases Brents algorithm needs slightly less and in other cases slightly more function evaluations).

• Julia package Roots.jl provides various algorithms to solve the problem above but it seems that Brents algorithm is not yet included (as of Nov. 2018).

```

module KINSOL - Solve nonlinear equation system with Sundials KINSOL

The goal is to solve the same system several times with KINSOL.