Source code for naginterfaces.library.examples.roots.sys_func_rcomm_ex

#!/usr/bin/env python3
"``naginterfaces.library.roots.sys_func_rcomm`` Python Example."

# NAG Copyright 2017-2019.

# pylint: disable=invalid-name,too-many-locals

from math import sqrt
import warnings

import numpy as np

from naginterfaces.base import utils
from naginterfaces.library import machine, roots

def main():
    """
    Example for :func:`naginterfaces.library.roots.sys_func_rcomm`.

    Find a solution of a system of nonlinear equations.

    Demonstrates promoting a ``NagAlgorithmicWarning`` to an
    exception.

    >>> main()
    naginterfaces.library.roots.sys_func_rcomm Python Example Results.
    Solution of a system of nonlinear equations.
    Final 2-norm of the residuals after 11 iterations is 1.19264e-08.
    Final approximate solution is
    (
      -0.5707,
      -0.6816,
      -0.7017,
      -0.7042,
      -0.7014,
      -0.6919,
      -0.6658,
      -0.5960,
      -0.4164,
    )
    """

    print('naginterfaces.library.roots.sys_func_rcomm Python Example Results.')
    print('Solution of a system of nonlinear equations.')

    # sys_func_rcomm categorizes three exit cases as 'Warnings',
    # for when progress to the solution stagnates.
    # warnings.simplefilter can be used to escalate these warnings to
    # exceptions:
    warnings.simplefilter('error', utils.NagAlgorithmicWarning)

    irevcm = 0

    # The initial value of x provides a rough solution.
    n = 9
    x = np.array([-1.]*n)
    xtol = sqrt(machine.precision())
    ml = 1
    mu = 1
    mode = 2
    diag = np.ones(n)

    # These arrays will be updated during the computation:
    fvec = np.empty(n)
    fjac = np.empty((n, n))
    r = np.empty(n*(n+1)//2)
    qtf = np.empty(n)

    # Communication object:
    comm = {}

    count = 0

    while True:
        irevcm = roots.sys_func_rcomm(
            irevcm, x, fvec, ml, mu, mode, diag, fjac, r, qtf, comm, xtol=xtol
        )

        if irevcm == 1:
            count += 1
            # Monitor progress here if desired.
            continue

        if irevcm == 2:
            # Evaluate the function.
            fvec = (3. - 2.*x)*x + 1.
            fvec[1:] = fvec[1:] - x[:n-1]
            fvec[:n-1] = fvec[:n-1] - 2.*x[1:]
            continue

        break

    norm_str = 'Final 2-norm of the residuals after {:d} iterations is {:.5e}.'
    print(norm_str.format(count, np.linalg.norm(fvec)))
    print('Final approximate solution is')
    print('(\n' + '\n'.join(['  {:.4f},']*len(x)).format(*x) + '\n)')

if __name__ == '__main__':
    import doctest
    import sys
    sys.exit(
        doctest.testmod(
            None, verbose=True, report=False,
            optionflags=doctest.REPORT_NDIFF,
        ).failed
    )