# Source code for naginterfaces.library.roots

# -*- coding: utf-8 -*-
r"""
Module Summary
--------------
Interfaces for the NAG Mark 29.2 roots Chapter.

roots - Roots of One or More Transcendental Equations

This module is concerned with the calculation of zeros of continuous functions of one or more variables.
The majority of problems considered are for real-valued functions of real variables, in which case complex equations must be expressed in terms of the equivalent larger system of real equations.

Background
----------
The module divides naturally into two parts.

A Single Equation
~~~~~~~~~~~~~~~~~
The first deals with the real zeros of a real function of a single variable :math:f\left(x\right).

There are three functions with simple calling sequences.
The first assumes that you can determine an initial interval :math:\left[a, b\right] within which the desired zero lies, (that is, where :math:f\left(a\right)\times f\left(b\right) < 0), and outside which all other zeros lie.
The function then systematically subdivides the interval to produce a final interval containing the zero.
This final interval has a length bounded by your specified error requirements; the end of the interval where the function has smallest magnitude is returned as the zero.
This function is guaranteed to converge to a **simple** zero of the function. (Here we define a simple zero as a zero corresponding to a sign-change of the function; none of the available functions are capable of making any finer distinction.) However, as with the other functions described below, a non-simple zero might be determined and it is left to you to check for this.
The algorithm used is due to Brent (1973).

The two other functions are both designed for the case where you are unable to specify an interval containing the simple zero.
One starts from an initial point and performs a search for an interval containing a simple zero.
If such an interval is computed then the method described above is used next to determine the zero accurately.
The other method uses a 'continuation' method based on a secant iteration.
A sequence of subproblems is solved; the first of these is trivial and the last is the actual problem of finding a zero of :math:f\left(x\right).
The intermediate problems employ the solutions of earlier problems to provide initial guesses for the secant iterations used to calculate their solutions.

Three other functions are also supplied.
They employ reverse communication and use the same core algorithms as the functions described above.

Finally, two functions are provided to return values of Lambert's :math:W function (sometimes known as the 'product log' or 'Omega' function), which is the inverse function of

.. math::

that is, if Lambert's :math:W function :math:W\left(x\right) = a for :math:\left. x, a\right. \in \mathbb{C}, then :math:a is a zero of the function :math:F\left(w\right) = we^w-x.
One function uses the iterative method described in Barry et al. (1995) to return values from the real branches of :math:W (restricting :math:\left. x, a\right. \in \mathbb{R}).
The second function enforces no such restriction, and uses the approach described in Corless et al. (1996).

Systems of Equations
~~~~~~~~~~~~~~~~~~~~
The functions in the second part of this module are designed to solve a set of nonlinear equations in :math:n unknowns

.. math::
f_i\left(x\right) = 0\text{, }\quad i = 1,2,\ldots,n\text{, }\quad x = {\left(x_1, x_2, \ldots, x_n\right)}^\mathrm{T}\text{,}

where :math:\mathrm{T} stands for transpose.

It is assumed that the functions are continuous and differentiable so that the matrix of first partial derivatives of the functions, the **Jacobian** matrix :math:J_{{ij}}\left(x\right) = \left(\frac{{\partial f_i}}{{\partial x_j}}\right) evaluated at the point :math:x, exists, though it may not be possible to calculate it directly.

The functions :math:f_i must be independent, otherwise there will be an infinity of solutions and the methods will fail.
However, even when the functions are independent the solutions may not be unique.
Since the methods are iterative, an initial guess at the solution has to be supplied, and the solution located will usually be the one closest to this initial guess.

Available Routines
------------------
Zeros of Functions of One Variable
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The functions can be divided into two classes.
There are three functions (:meth:contfn_interval_rcomm, :meth:contfn_cntin_rcomm and :meth:contfn_brent_rcomm) all written in reverse communication form and three (:meth:contfn_brent_interval, :meth:contfn_cntin and :meth:contfn_brent) written in direct communication form (see Direct and Reverse Communication Routines <https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/h/howtounl.html#revcommroutines>__ for a description of the difference between these two conventions).
The direct communication functions are designed for inexperienced users and, in particular, for solving problems where the function :math:f\left(x\right) whose zero is to be calculated, can be coded as a callback function.
These functions find the zero by using the same core algorithms as the reverse communication functions.
Experienced users are recommended to use the reverse communication functions directly as they permit you more control of the calculation.
Indeed, if the zero-finding process is embedded in a much larger program then the reverse communication functions should always be used.

The recommendation as to which function should be used depends mainly on whether you can supply an interval :math:\left[a, b\right] containing the zero; that is, where :math:f\left(a\right)\times f\left(b\right) < 0.
If the interval can be supplied, then :meth:contfn_brent (or, in reverse communication, :meth:contfn_brent_rcomm) should be used, in general.
This recommendation should be qualified in the case when the only interval which can be supplied is very long relative to your error requirements **and** you can also supply a good approximation to the zero.
In this case :meth:contfn_cntin (or, in reverse communication, :meth:contfn_cntin_rcomm) **may** prove more efficient (though these latter functions will not provide the error bound available from :meth:contfn_brent_rcomm).

If an interval containing the zero cannot be supplied then you must choose between :meth:contfn_brent_interval (or, in reverse communication, :meth:contfn_interval_rcomm followed by :meth:contfn_brent_rcomm) and :meth:contfn_cntin (or, in reverse communication, :meth:contfn_cntin_rcomm). :meth:contfn_brent_interval first determines an interval containing the zero, and then proceeds as in :meth:contfn_brent; it is particularly recommended when you do not have a good initial approximation to the zero.
If a good initial approximation to the zero is available then :meth:contfn_cntin is to be preferred.
Since neither of these latter functions has guaranteed convergence to the zero, you are recommended to experiment with both in case of difficulty.

Solution of Sets of Nonlinear Equations
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The solution of a set of nonlinear equations

.. math::
f_i\left(x_1, x_2, \ldots, x_n\right) = 0\text{, }\quad i = 1,2,\ldots,n

can be regarded as a special case of the problem of finding a minimum of a sum of squares

.. math::
s\left(x\right) = \sum_{{i = 1}}^m\left[f_i\left(x_1, x_2, \ldots, x_n\right)\right]^2\text{, }\quad \left(m\geq n\right)\text{.}

So the functions in submodule :mod:~naginterfaces.library.opt are relevant as well as the special nonlinear equations functions.

The functions for solving a set of nonlinear equations can also be divided into classes.
There are six functions (:meth:sys_func_aa, :meth:sys_func_easy, :meth:sys_func_expert, :meth:sparsys_func_easy, :meth:sys_deriv_easy and :meth:sys_deriv_expert) all written in direct communication form and three (:meth:sys_func_aa_rcomm, :meth:sys_func_rcomm and :meth:sys_deriv_rcomm) written in reverse communication form.
The direct communication functions are designed for inexperienced users and, in particular, these functions require the :math:f_i (and possibly their derivatives) to be calculated in functions.
These should be set up carefully so the NAG Library functions can work as efficiently as possible.
Experienced users are recommended to use the reverse communication functions as they permit you more control of the calculation.
Indeed, if the zero-finding process is embedded in a much larger program then the reverse communication functions should always be used.

The main decision you have to make is whether to supply the derivatives :math:\frac{{\partial f_i}}{{\partial x_j}}.
It is advisable to do so if possible, since the results obtained by algorithms which use derivatives are generally more reliable than those obtained by algorithms which do not use derivatives.

:meth:sys_deriv_easy, :meth:sys_deriv_expert and :meth:sys_deriv_rcomm require you to provide the derivatives, whilst :meth:sys_func_aa, :meth:sys_func_aa_rcomm, :meth:sys_func_easy, :meth:sys_func_expert, :meth:sys_func_rcomm and :meth:sparsys_func_easy do not. :meth:sys_func_easy, :meth:sys_deriv_easy and :meth:sparsys_func_easy are easy-to-use functions; greater flexibility may be obtained using :meth:sys_func_expert and :meth:sys_deriv_expert (or, in reverse communication, :meth:sys_func_rcomm and :meth:sys_deriv_rcomm), but these have longer argument lists. :meth:sys_func_easy, :meth:sys_func_expert, :meth:sys_func_rcomm and :meth:sparsys_func_easy approximate the derivatives internally using finite differences. :meth:sys_func_aa and :meth:sys_func_aa_rcomm do not use derivatives at all, and may be useful when the cost of evaluating :math:f\left(x\right) is high.
If :math:f\left(x\right) can be evaluated cheaply, then functions which use the Jacobian or its approximation may converge more quickly.

:meth:sparsys_func_easy is an easy-to-use function specially adapted for sparse problems, that is, problems where each function depends on a small subset of the :math:n variables so that the Jacobian matrix has many zeros.
It employs sparse linear algebra methods and consequently is expected to take significantly less time to complete than the other functions, especially if :math:n is large.

:meth:sys_deriv_check is provided for use in conjunction with :meth:sys_deriv_easy, :meth:sys_deriv_expert and :meth:sys_deriv_rcomm to check the user-supplied derivatives for consistency with the functions themselves.
You are strongly advised to make use of this function whenever :meth:sys_deriv_easy, :meth:sys_deriv_expert or :meth:sys_deriv_rcomm is used.

Firstly, the calculation of the functions and their derivatives should be ordered so that **cancellation errors** are avoided.
This is particularly important in a function that uses these quantities to build up estimates of higher derivatives.

Secondly, **scaling** of the variables has a considerable effect on the efficiency of a function.
The problem should be designed so that the elements of :math:x are of similar magnitude.
The same comment applies to the functions, i.e., all the :math:f_i should be of comparable size.

The accuracy is usually determined by the accuracy arguments of the functions, but the following points may be useful.

(i) Greater accuracy in the solution may be requested by choosing smaller input values for the accuracy arguments. However, if unreasonable accuracy is demanded, rounding errors may become important and cause a failure.

(#) Some idea of the accuracies of the :math:x_i may be obtained by monitoring the progress of the function to see how many figures remain unchanged during the last few iterations.

(#) An approximation to the error in the solution :math:x is given by :math:e where :math:e is the solution to the set of linear equations

.. math::
J\left(x\right)e = -f\left(x\right)

where :math:f\left(x\right) = {\left({f_1\left(x\right)}, {f_2\left(x\right)}, \ldots, {f_n\left(x\right)}\right)}^\mathrm{T}.

Note that the :math:QR decomposition of :math:J is available from :meth:sys_func_expert and :meth:sys_deriv_expert (or, in reverse communication, :meth:sys_func_rcomm and :meth:sys_deriv_rcomm) so that

.. math::
Re = -Q^\mathrm{T}f

and :math:Q^\mathrm{T}f is also provided by these functions.

(#) If the functions :math:f_i\left(x\right) are changed by small amounts :math:\epsilon_i, for :math:\textit{i} = 1,2,\ldots,n, then the corresponding change in the solution :math:x is given approximately by :math:\sigma, where :math:\sigma is the solution of the set of linear equations

.. math::
J\left(x\right)\sigma = -\epsilon \text{.}

Thus one can estimate the sensitivity of :math:x to any uncertainties in the specification of :math:f_i\left(x\right), for :math:\textit{i} = 1,2,\ldots,n.
As noted above, the sophisticated functions :meth:sys_func_expert and :meth:sys_deriv_expert (or, in reverse communication, :meth:sys_func_rcomm and :meth:sys_deriv_rcomm) provide the :math:QR decomposition of :math:J.

Values of Lambert's :math:W function
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
If you require purely-real values of :math:W, these will be evaluated marginally more efficiently by :meth:lambertw_real than :meth:lambertw_complex owing to the differing iterative procedures used by each function.

Tutorial
--------
Zeros of Functions of One Variable
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Consider the problem of finding :math:x such that :math:x = \mathrm{exp}\left({-x}\right); i.e., the value of :math:x that gives :math:f(x) = 0 when :math:f(x) = x-\mathrm{exp}\left({-x}\right).

>>> from math import exp
>>> f = lambda x: x - exp(-x)

All of the NAG zero-finding routines require you to have some basic knowledge of the structure of the function you are addressing in order that the solver can know where to begin.
This knowledge must either be in the form of a rough initial guess for the zero, or as an interval on which :math:f changes sign (and, therefore, passes through zero).

In the worst case a sign-change interval is not known and the initial guess for the zero may be fairly poor:

>>> x = 1.

If so, use :meth:contfn_brent_interval:

>>> from naginterfaces.library import roots
>>> soln = roots.contfn_brent_interval(x, f)

Upon successful exit :math:\texttt{soln} contains the computed zero and the bounding interval that was used internally:

>>> print('Final value of x is {:1.8f}.'.format(soln.x))
Final value of x is 0.56714329.
>>> print(
...     'Search interval was [' +
...     ', '.join(['{:1.8f}']*2).format(soln.a, soln.b) +
...     '].'
... )
Search interval was [0.50000000, 0.90000000].

For an alternative approach that uses continuation, more suitable for when a better initial guess is known, use :meth:contfn_cntin:

>>> x = 0.6
>>> soln = roots.contfn_cntin(x, f)

The continuation algorithm returns only the final :math:x:

>>> print('Final value of x is {:1.8f}.'.format(soln))
Final value of x is 0.56714329.

Solution of Sets of Nonlinear Equations
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Consider the following tridiagonal system of nine equations in nine variables:

.. math::
\begin{array}{rcl}\left(3-2x_1\right)x_1-2x_2& = &-1\text{,}\\-x_{{i-1}}+\left(3-2x_i\right)x_i-2x_{{i+1}}& = &-1\text{, }\quad i = 2,3,\ldots,8,\\-x_8+\left(3-2x_9\right)x_9& = &-1\text{.}\end{array}

as a vector function

>>> def f(x):
...     fvec = (3. - 2.*x)*x + 1.
...     fvec[1:] = fvec[1:] - x[:n-1]
...     fvec[:n-1] = fvec[:n-1] - 2.*x[1:]
...     return fvec

for which we seek :math:x such that :math:f(x) = 0.

As with scalar root-finding above, the NAG solver needs you to supply a rough initial guess for :math:x.

>>> import numpy as np
>>> n = 9; x = -1.*np.ones(n)

To solve a system of nonlinear equations, again the NAG routines can benefit from your knowledge of the structure of the problem as outlined above.

If you are happy for the Jacobian of the system to be approximated internally, use, say, :meth:sys_func_easy:

>>> from naginterfaces.library import roots
>>> soln = roots.sys_func_easy(f, x)
>>> def show_soln():
...     print(
...         'Final approximate solution is\n'
...         '(\n' +
...         '\n'.join(['  {:.4f},']*n).format(*soln.x) +
...         '\n)'
...     )
>>> show_soln()
Final approximate solution is
(
-0.5707,
-0.6816,
-0.7017,
-0.7042,
-0.7014,
-0.6919,
-0.6658,
-0.5960,
-0.4164,
)

The exact Jacobian of the system in question can be represented using the following Python function:

>>> import numpy as np
>>> def fjac(x):
...     fj = np.zeros((len(x), len(x)))
...     fj[0, 0] = 3. - 4.*x
...     fj[0, 1] = -2.
...     for i in range(1, n-2):
...         fj[i, i-1] = -1.
...         fj[i, i] = 3. - 4.*x[i]
...         fj[i, i+1] = -2.
...     fj[n-1, n-2] = -1.
...     fj[n-1, n-1] = 3. - 4.*x[n-1]
...     return fj

Function :meth:sys_deriv_easy can be used to take advantage of these exact derivatives:

>>> from naginterfaces.library import roots
>>> soln = roots.sys_deriv_easy(f, fjac, x)
>>> show_soln()
Final approximate solution is
(
-0.5707,
-0.6816,
-0.7017,
-0.7042,
-0.7014,
-0.6919,
-0.6658,
-0.5960,
-0.4164,
)

When the Jacobian is known to contain a large number of zero entries, but isn't available in an exact form, :meth:sparsys_func_easy may be appropriate.
With this solver we are asked to compute values of the system for a supplied base-:math:1 index mask. (In the contrived example below we use precomputed 'full' values for the system, using our definition of :math:\texttt{f} above, to populate the requested entries of the return vector.
For a real-world sparse problem where memory usage is a concern this would clearly be counterproductive.) The function vector is updated in place in the callback, which further optimizes memory usage.

>>> def sparsf(indf, x, fvec):
...     fv = f(x)
...     for ind in indf:
...        fvec[ind-1] = fv[ind-1]

For our system the number of nonzeros is :math:n+2\left(n-1\right).
Communicating this in the NAG call can reduce memory requirements.

>>> from naginterfaces.library import roots
>>> soln = roots.sparsys_func_easy(sparsf, x, nnz=3*n-2)
>>> show_soln()
Final approximate solution is
(
-0.5707,
-0.6816,
-0.7017,
-0.7042,
-0.7014,
-0.6919,
-0.6658,
-0.5960,
-0.4164,
)

Computational efficiency may also be enhanced by using the fixed-point iteration routine :meth:sys_func_aa and its Anderson acceleration facility.
The potential downsides are that convergence is less reliable than with the Powell hybrid method used in the solvers shown previously, and that the system may need to be rearranged into the form :math:g(x)-x = 0 for best performance:

>>> import numpy as np
>>> def fixedf(x):
...     fvec = np.zeros(n)
...     fvec = (2.*x - 1.)/(3. - 2.*x)
...     for i in range(1, n-1):
...         fvec[i] = (2.*x[i+1] - 1. + x[i-1])/(3. - 2.*x[i])
...     fvec[n-1] = (x[n-2] - 1.)/(3. - 2.*x[n-1])
...     return fvec - x

The size of the iteration subspace cache to use is supplied via argument :math:\texttt{m}:

>>> from naginterfaces.library import roots
>>> soln = roots.sys_func_aa(fixedf, x, m=4)
>>> show_soln()
Final approximate solution is
(
-0.5707,
-0.6816,
-0.7017,
-0.7042,
-0.7014,
-0.6919,
-0.6658,
-0.5960,
-0.4164,
)

--------
naginterfaces.library.examples.roots :
This subpackage contains examples for the roots module.
See also the :ref:library_roots_ex subsection.

Functionality Index
-------------------

**Lambert's** :math:W **function**

complex values: :meth:lambertw_complex

real values: :meth:lambertw_real

**Zeros of functions of one variable**

direct communication

binary search followed by Brent algorithm: :meth:contfn_brent_interval

Brent algorithm: :meth:contfn_brent

continuation method: :meth:contfn_cntin

reverse communication

binary search: :meth:contfn_interval_rcomm

Brent algorithm: :meth:contfn_brent_rcomm

continuation method: :meth:contfn_cntin_rcomm

**Zeros of functions of several variables**

checking function

checks user-supplied Jacobian: :meth:sys_deriv_check

direct communication

Anderson acceleration: :meth:sys_func_aa

easy-to-use

derivatives required: :meth:sys_deriv_easy

no derivatives required: :meth:sys_func_easy

no derivatives required, sparse: :meth:sparsys_func_easy

sophisticated

derivatives required: :meth:sys_deriv_expert

no derivatives required: :meth:sys_func_expert

reverse communication

Anderson acceleration: :meth:sys_func_aa_rcomm

sophisticated

derivatives required: :meth:sys_deriv_rcomm

no derivatives required: :meth:sys_func_rcomm

For full information please refer to the NAG Library document

https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05intro.html
"""

from __future__ import absolute_import

from math import sqrt

from ..library import machine

[docs]def contfn_brent_interval(x, f, h=0.1, eps=1000.0*machine.precision(), eta=0.0, data=None):
r"""
contfn_brent_interval locates a simple zero of a continuous function from a given starting value.
It uses a binary search to locate an interval containing a zero of the function, then Brent's method, which is a combination of nonlinear interpolation, linear extrapolation and bisection, to locate the zero precisely.

.. _c05au-py2-py-doc:

For full information please refer to the NAG Library document for c05au

https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05auf.html

.. _c05au-py2-py-parameters:

**Parameters**
**x** : float
An initial approximation to the zero.

**f** : callable retval = f(x, data=None)
:math:\mathrm{f} must evaluate the function :math:f whose zero is to be determined.

**Parameters**
**x** : float
The point at which the function must be evaluated.

**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.

**Returns**
**retval** : float
The value of :math:f at the specified point.

**h** : float, optional
A step length for use in the binary search for an interval containing the zero. The maximum interval searched is :math:\left[{\mathrm{x}-256.0\times \mathrm{h}}, {\mathrm{x}+256.0\times \mathrm{h}}\right].

**eps** : float, optional
The termination tolerance on :math:x (see :ref:Notes <c05au-py2-py-notes>).

**eta** : float, optional
A value such that if :math:\left\lvert f\left(x\right)\right\rvert \leq \mathrm{eta}, :math:x is accepted as the zero. :math:\mathrm{eta} may be specified as :math:0.0 (see Accuracy <https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05auf.html#accuracy>__).

**data** : arbitrary, optional
User-communication data for callback functions.

**Returns**
**x** : float
If the function exits successfully or :math:\mathrm{errno} = 4, :math:\mathrm{x} is the final approximation to the zero.

If :math:\mathrm{errno} = 3, :math:\mathrm{x} is likely to be a pole of :math:f\left(x\right).

Otherwise, :math:\mathrm{x} contains no useful information.

**a** : float
The lower bound of the interval resulting from the binary search

**b** : float
The upper bound of the interval resulting from the binary search

.. _c05au-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:1)
On entry, :math:\mathrm{eps} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{eps} > 0.0.

(errno :math:1)
On entry, :math:\mathrm{x} = \langle\mathit{\boldsymbol{value}}\rangle and :math:\mathrm{h} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:{\mathrm{x}+\mathrm{h}} \neq \mathrm{x} (to machine accuracy).

(errno :math:2)
An interval containing the zero could not be found. Increasing :math:\mathrm{h} and calling contfn_brent_interval again will increase the range searched for the zero. Decreasing :math:\mathrm{h} and calling contfn_brent_interval again will refine the mesh used in the search for the zero.

**Warns**
**NagAlgorithmicWarning**
(errno :math:3)
Solution may be a pole rather than a zero.

(errno :math:4)
The tolerance :math:\mathrm{eps} has been set too small for the problem being solved. However, the value :math:\mathrm{x} returned is a good approximation to the zero. :math:\mathrm{eps} = \langle\mathit{\boldsymbol{value}}\rangle.

.. _c05au-py2-py-notes:

**Notes**
contfn_brent_interval attempts to locate an interval :math:\left[a, b\right] containing a simple zero of the function :math:f\left(x\right) by a binary search starting from the initial point :math:x = \mathrm{x} and using repeated calls to :meth:contfn_interval_rcomm.
If this search succeeds, then the zero is determined to a user-specified accuracy by a call to :meth:contfn_brent.
The specifications of functions :meth:contfn_interval_rcomm and :meth:contfn_brent should be consulted for details of the methods used.

The approximation :math:x to the zero :math:\alpha is determined so that at least one of the following criteria is satisfied:

(i) :math:\left\lvert x-\alpha \right\rvert \leq \mathrm{eps},

(#) :math:\left\lvert f\left(x\right)\right\rvert \leq \mathrm{eta}.

.. _c05au-py2-py-references:

**References**
Brent, R P, 1973, Algorithms for Minimization Without Derivatives, Prentice--Hall

--------
:meth:naginterfaces.library.examples.roots.contfn_brent_interval_ex.main
"""
raise NotImplementedError

[docs]def contfn_interval_rcomm(x, fx, h, boundl, boundu, y, comm, ind):
r"""
contfn_interval_rcomm attempts to locate an interval containing a simple zero of a continuous function using a binary search.
It uses reverse communication for evaluating the function.

.. _c05av-py2-py-doc:

For full information please refer to the NAG Library document for c05av

https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05avf.html

.. _c05av-py2-py-parameters:

**Parameters**
**x** : float
On initial entry: the best available approximation to the zero.

**fx** : float
On initial entry: if :math:\mathrm{ind} = 1, :math:\mathrm{fx} need not be set.

If :math:\mathrm{ind} = -1, :math:\mathrm{fx} must contain :math:f\left(\mathrm{x}\right) for the initial value of :math:\mathrm{x}.

On intermediate entry: must contain :math:f\left(\mathrm{x}\right) for the current value of :math:\mathrm{x}.

**h** : float
On initial entry: a basic step size which is used in the binary search for an interval containing a zero. The basic step sizes :math:\mathrm{h},0.1\times \mathrm{h}, :math:0.01\times \mathrm{h} and :math:0.001\times \mathrm{h} are used in turn when searching for the zero.

**boundl** : float
The lower bound for the interval of search for the zero.

**boundu** : float
The upper bound for the interval of search for the zero.

**y** : float
On initial entry: need not be set.

**comm** : dict, communication object, modified in place
Communication structure.

On initial entry: need not be set.

**ind** : int
On initial entry: must be set to :math:1 or :math:-1.

:math:\mathrm{ind} = 1

:math:\mathrm{fx} need not be set.

:math:\mathrm{ind} = -1

:math:\mathrm{fx} must contain :math:f\left(\mathrm{x}\right).

**Returns**
**x** : float
On intermediate exit: contains the point at which :math:f must be evaluated before re-entry to the function.

On final exit: contains one end of an interval containing the zero, the other end being in :math:\mathrm{y}, unless an error has occurred. If :math:\mathrm{errno} = 4, :math:\mathrm{x} and :math:\mathrm{y} are the end points of the largest interval searched. If a zero is located exactly, its value is returned in :math:\mathrm{x} (and in :math:\mathrm{y}).

**h** : float
On final exit: is undefined.

**y** : float
On final exit: contains the closest point found to the final value of :math:\mathrm{x}, such that :math:f\left(\mathrm{x}\right)\times f\left(\mathrm{y}\right)\leq 0.0. If a value :math:\mathrm{x} is found such that :math:f\left(\mathrm{x}\right) = 0, :math:\mathrm{y} = \mathrm{x}. On final exit with :math:\mathrm{errno} = 4, :math:\mathrm{x} and :math:\mathrm{y} are the end points of the largest interval searched.

**c1** : float
On final exit: if the function exits successfully or :math:\mathrm{errno} = 4, :math:\mathrm{c1} is set to :math:f\left(\mathrm{y}\right).

**ind** : int
On intermediate exit: contains :math:2 or :math:3. The calling program must evaluate :math:f at :math:\mathrm{x}, storing the result in :math:\mathrm{fx}, and re-enter contfn_interval_rcomm with all other arguments unchanged.

On final exit: contains :math:0.

.. _c05av-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:1)
On entry, :math:\mathrm{boundl} = \langle\mathit{\boldsymbol{value}}\rangle and :math:\mathrm{boundu} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{boundl} < \mathrm{boundu}.

(errno :math:1)
On entry, :math:\mathrm{x} = \langle\mathit{\boldsymbol{value}}\rangle, :math:\mathrm{boundl} = \langle\mathit{\boldsymbol{value}}\rangle and :math:\mathrm{boundu} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{boundl}\leq \mathrm{x}\leq \mathrm{boundu}.

(errno :math:1)
On entry, :math:\mathrm{x} = \langle\mathit{\boldsymbol{value}}\rangle, :math:\mathrm{h} = \langle\mathit{\boldsymbol{value}}\rangle, :math:\mathrm{boundl} = \langle\mathit{\boldsymbol{value}}\rangle and :math:\mathrm{boundu} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: either :math:\mathrm{x}+\mathrm{h} or :math:\mathrm{x}-\mathrm{h} must lie inside the closed interval :math:\left[\mathrm{boundl}, \mathrm{boundu}\right].

(errno :math:2)
On entry, :math:\mathrm{x} = \langle\mathit{\boldsymbol{value}}\rangle and :math:\mathrm{h} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{h} must be sufficiently large that :math:{\mathrm{x}+\mathrm{h}} \neq \mathrm{x} on the computer.

(errno :math:3)
On entry, :math:\mathrm{ind} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{ind} = -1, :math:1, :math:2 or :math:3.

**Warns**
**NagAlgorithmicWarning**
(errno :math:4)
An interval containing the zero could not be found. Try restarting with modified :math:\mathrm{x} and :math:\mathrm{h}.

.. _c05av-py2-py-notes:

**Notes**
You must supply an initial point :math:\mathrm{x} and a step :math:\mathrm{h}. contfn_interval_rcomm attempts to locate a short interval :math:\left[\mathrm{x}, \mathrm{y}\right]⊂\left[\mathrm{boundl}, \mathrm{boundu}\right] containing a simple zero of :math:f\left(x\right).

(On exit we may have :math:\mathrm{x} > \mathrm{y}; :math:\mathrm{x} is determined as the first point encountered in a binary search where the sign of :math:f\left(x\right) differs from the sign of :math:f\left(x\right) at the initial input point :math:\mathrm{x}.) The function attempts to locate a zero of :math:f\left(x\right) using :math:\mathrm{h}, :math:0.1\times \mathrm{h}, :math:0.01\times \mathrm{h} and :math:0.001\times \mathrm{h} in turn as its basic step before quitting with an error exit if unsuccessful.

contfn_interval_rcomm returns to the calling program for each evaluation of :math:f\left(x\right).
On each return you should set :math:\mathrm{fx} = f\left(\mathrm{x}\right) and call contfn_interval_rcomm again.
"""
raise NotImplementedError

[docs]def contfn_cntin(x, f, eps=1000.0*machine.precision(), eta=0.0, nfmax=1000, data=None):
r"""
contfn_cntin attempts to locate a zero of a continuous function using a continuation method based on a secant iteration.

.. _c05aw-py2-py-doc:

For full information please refer to the NAG Library document for c05aw

https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05awf.html

.. _c05aw-py2-py-parameters:

**Parameters**
**x** : float
An initial approximation to the zero.

**f** : callable retval = f(x, data=None)
:math:\mathrm{f} must evaluate the function :math:f whose zero is to be determined.

**Parameters**
**x** : float
The point at which the function must be evaluated.

**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.

**Returns**
**retval** : float
The value of :math:f evaluated at :math:\mathrm{x}.

**eps** : float, optional
An absolute tolerance to control the accuracy to which the zero is determined. In general, the smaller the value of :math:\mathrm{eps} the more accurate :math:\mathrm{x} will be as an approximation to :math:\alpha. Indeed, for very small positive values of :math:\mathrm{eps}, it is likely that the final approximation will satisfy :math:\left\lvert \mathrm{x}-\alpha \right\rvert < \mathrm{eps}. You are advised to call the function with more than one value for :math:\mathrm{eps} to check the accuracy obtained.

**eta** : float, optional
A value such that if :math:\left\lvert f\left(x\right)\right\rvert < \mathrm{eta}, :math:x is accepted as the zero. :math:\mathrm{eta} may be specified as :math:0.0 (see Accuracy <https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05awf.html#accuracy>__).

**nfmax** : int, optional
The maximum permitted number of calls to :math:\mathrm{f} from contfn_cntin. If :math:\mathrm{f} is inexpensive to evaluate, :math:\mathrm{nfmax} should be given a large value (say :math:\text{} > 1000).

**data** : arbitrary, optional
User-communication data for callback functions.

**Returns**
**x** : float
The final approximation to the zero, unless :math:\mathrm{errno} = 1 or 2, in which case it contains no useful information.

.. _c05aw-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:1)
On entry, :math:\mathrm{nfmax} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{nfmax} > 0.

(errno :math:1)
On entry, :math:\mathrm{eps} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{eps} > 0.0.

(errno :math:2)
Internal scale factor invalid for this problem. Consider using :meth:contfn_cntin_rcomm instead and setting :math:\textit{scal}.

(errno :math:3)
Either :math:\mathrm{f} has no zero near :math:\mathrm{x} or too much accuracy has been requested. Check the coding of :math:\mathrm{f} or increase :math:\mathrm{eps}.

(errno :math:4)
More than :math:\mathrm{nfmax} calls have been made to :math:\mathrm{f}.

.. _c05aw-py2-py-notes:

**Notes**
contfn_cntin attempts to obtain an approximation to a simple zero :math:\alpha of the function :math:f\left(x\right) given an initial approximation :math:x to :math:\alpha.
The zero is found by a call to :meth:contfn_cntin_rcomm whose specification should be consulted for details of the method used.

The approximation :math:x to the zero :math:\alpha is determined so that at least one of the following criteria is satisfied:

(i) :math:\left\lvert x-\alpha \right\rvert \sim \mathrm{eps},

(#) :math:\left\lvert f\left(x\right)\right\rvert < \mathrm{eta}.
"""
raise NotImplementedError

[docs]def contfn_cntin_rcomm(x, fx, tol, ir, comm, ind, scal=sqrt(machine.precision())):
r"""
contfn_cntin_rcomm attempts to locate a zero of a continuous function using a continuation method based on a secant iteration.
It uses reverse communication for evaluating the function.

.. _c05ax-py2-py-doc:

For full information please refer to the NAG Library document for c05ax

https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05axf.html

.. _c05ax-py2-py-parameters:

**Parameters**
**x** : float
On initial entry: an initial approximation to the zero.

**fx** : float
On initial entry: if :math:\mathrm{ind} = 1, :math:\mathrm{fx} need not be set.

If :math:\mathrm{ind} = -1, :math:\mathrm{fx} must contain :math:f\left(\mathrm{x}\right) for the initial value of :math:\mathrm{x}.

On intermediate entry: must contain :math:f\left(\mathrm{x}\right) for the current value of :math:\mathrm{x}.

**tol** : float
On initial entry: a value that controls the accuracy to which the zero is determined. :math:\mathrm{tol} is used in determining the convergence of the secant iteration used at each stage of the continuation process. It is used directly when solving the last problem (:math:\theta_m = 0 in :ref:Notes <c05ax-py2-py-notes>), and :math:\sqrt{\mathrm{tol}} is used for the problem defined by :math:\theta_r, :math:r < m. Convergence to the accuracy specified by :math:\mathrm{tol} is not guaranteed, and so you are recommended to find the zero using at least two values for :math:\mathrm{tol} to check the accuracy obtained.

**ir** : int
On initial entry: indicates the type of error test required, as follows. Solving the problem defined by :math:\theta_r, :math:1\leq r\leq m, involves computing a sequence of secant iterates :math:x_r^0,x_r^1,\ldots \text{}. This sequence will be considered to have converged only if:

for :math:\mathrm{ir} = 0,

.. math::
\left\lvert x_r^{\left(i+1\right)}-x_r^{\left(i\right)}\right\rvert \leq \textit{eps}\times \mathrm{max}\left(1.0, \left\lvert x_r^{\left(i\right)}\right\rvert \right)\text{,}

for :math:\mathrm{ir} = 1,

.. math::
\left\lvert x_r^{\left(i+1\right)}-x_r^{\left(i\right)}\right\rvert \leq \textit{eps}\text{,}

for :math:\mathrm{ir} = 2,

.. math::
\left\lvert x_r^{\left(i+1\right)}-x_r^{\left(i\right)}\right\rvert \leq \textit{eps}\times \left\lvert x_r^{\left(i\right)}\right\rvert \text{,}

for some :math:i > 1; here :math:\textit{eps} is either :math:\mathrm{tol} or :math:\sqrt{\mathrm{tol}} as discussed above.

Note that there are other subsidiary conditions (not given here) which must also be satisfied before the secant iteration is considered to have converged.

**comm** : dict, communication object, modified in place
Communication structure.

On initial entry: need not be set.

**ind** : int
On initial entry: must be set to :math:1 or :math:-1.

:math:\mathrm{ind} = 1

:math:\mathrm{fx} need not be set.

:math:\mathrm{ind} = -1

:math:\mathrm{fx} must contain :math:f\left(\mathrm{x}\right).

**scal** : float, optional
On initial entry: a factor for use in determining a significant approximation to the derivative of :math:f\left(x\right) at :math:x = x_0, the initial value. A number of difference approximations to :math:f^{\prime }\left(x_0\right) are calculated using

.. math::
f^{\prime }\left(x_0\right)\sim \left(f\left(x_0+h\right)-f\left(x_0\right)\right)/h

where :math:\left\lvert h\right\rvert < \left\lvert \mathrm{scal}\right\rvert and :math:h has the same sign as :math:\mathrm{scal}. A significance (cancellation) check is made on each difference approximation and the approximation is rejected if insignificant.

**Returns**
**x** : float
On intermediate exit: the point at which :math:f must be evaluated before re-entry to the function.

On final exit: the final approximation to the zero.

**ind** : int
On intermediate exit: contains :math:2, :math:3 or :math:4. The calling program must evaluate :math:f at :math:\mathrm{x}, storing the result in :math:\mathrm{fx}, and re-enter contfn_cntin_rcomm with all other arguments unchanged.

On final exit: contains :math:0.

.. _c05ax-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:1)
On entry, :math:\mathrm{tol} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{tol} > 0.0.

(errno :math:1)
On entry, :math:\mathrm{ir} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{ir} = 0, :math:1 or :math:2.

(errno :math:2)
On initial entry, :math:\mathrm{ind} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{ind} = -1 or :math:1.

(errno :math:2)
On intermediate entry, :math:\mathrm{ind} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{ind} = 2, :math:3 or :math:4.

(errno :math:3)
On entry, :math:\mathrm{scal} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:{\mathrm{x}+\mathrm{scal}} \neq \mathrm{x} (to machine accuracy).

(errno :math:3)
Significant derivatives of :math:f cannot be computed. This can happen when :math:f is almost constant and nonzero, for any value of :math:\mathrm{scal}.

(errno :math:4)
Current problem in the continuation sequence cannot be solved. Perhaps the original problem had no solution or the continuation path passes through a set of insoluble problems: consider refining the initial approximation to the zero. Alternatively, :math:\mathrm{tol} is too small, and the accuracy requirement is too stringent, or too large and the initial approximation too poor.

(errno :math:5)
Continuation away from the initial point is not possible. This error exit will usually occur if the problem has not been properly posed or the error requirement is extremely stringent.

(errno :math:6)
Final problem (with :math:\theta_m = 0) cannot be solved. It is likely that too much accuracy has been requested, or that the zero is at :math:\alpha = 0 and :math:\mathrm{ir} = 2.

.. _c05ax-py2-py-notes:

**Notes**
contfn_cntin_rcomm uses a modified version of an algorithm given in Swift and Lindfield (1978) to compute a zero :math:\alpha of a continuous function :math:f\left(x\right).
The algorithm used is based on a continuation method in which a sequence of problems

.. math::

are solved, where :math:1 = \theta_0 > \theta_1 > \cdots > \theta_m = 0 (the value of :math:m is determined as the algorithm proceeds) and where :math:x_0 is your initial estimate for the zero of :math:f\left(x\right).
For each :math:\theta_r the current problem is solved by a robust secant iteration using the solution from earlier problems to compute an initial estimate.

You must supply an error tolerance :math:\mathrm{tol}. :math:\mathrm{tol} is used directly to control the accuracy of solution of the final problem (:math:\theta_m = 0) in the continuation method, and :math:\sqrt{\mathrm{tol}} is used to control the accuracy in the intermediate problems (:math:\theta_1,\theta_2,\ldots,\theta_{{m-1}}).

.. _c05ax-py2-py-references:

**References**
Swift, A and Lindfield, G R, 1978, Comparison of a continuation method for the numerical solution of a single nonlinear equation, Comput. J. (21), 359--362
"""
raise NotImplementedError

[docs]def contfn_brent(a, b, f, eps=1000.0*machine.precision(), eta=0.0, data=None):
r"""
contfn_brent locates a simple zero of a continuous function in a given interval using Brent's method, which is a combination of nonlinear interpolation, linear extrapolation and bisection.

.. _c05ay-py2-py-doc:

For full information please refer to the NAG Library document for c05ay

https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05ayf.html

.. _c05ay-py2-py-parameters:

**Parameters**
**a** : float
:math:a, the lower bound of the interval.

**b** : float
:math:b, the upper bound of the interval.

**f** : callable retval = f(x, data=None)
:math:\mathrm{f} must evaluate the function :math:f whose zero is to be determined.

**Parameters**
**x** : float
The point at which the function must be evaluated.

**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.

**Returns**
**retval** : float
The value of :math:f evaluated at :math:\mathrm{x}.

**eps** : float, optional
The termination tolerance on :math:x (see :ref:Notes <c05ay-py2-py-notes>).

**eta** : float, optional
A value such that if :math:\left\lvert f\left(x\right)\right\rvert \leq \mathrm{eta}, :math:x is accepted as the zero. :math:\mathrm{eta} may be specified as :math:0.0 (see Accuracy <https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05ayf.html#accuracy>__).

**data** : arbitrary, optional
User-communication data for callback functions.

**Returns**
**x** : float
If the function exits successfully or :math:\mathrm{errno} = 2, :math:\mathrm{x} is the final approximation to the zero. If :math:\mathrm{errno} = 3, :math:\mathrm{x} is likely to be a pole of :math:f\left(x\right). Otherwise, :math:\mathrm{x} contains no useful information.

.. _c05ay-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:1)
On entry, :math:\mathrm{eps} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{eps} > 0.0.

(errno :math:1)
On entry, :math:\mathrm{a} = \langle\mathit{\boldsymbol{value}}\rangle and :math:\mathrm{b} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{a} \neq \mathrm{b}.

(errno :math:1)
On entry, :math:\mathrm{f}\left(\mathrm{a}\right) and :math:\mathrm{f}\left(\mathrm{b}\right) have the same sign with neither equalling :math:0.0: :math:\mathrm{f}\left(\mathrm{a}\right) = \langle\mathit{\boldsymbol{value}}\rangle and :math:\mathrm{f}\left(\mathrm{b}\right) = \langle\mathit{\boldsymbol{value}}\rangle.

**Warns**
**NagAlgorithmicWarning**
(errno :math:2)
No further improvement in the solution is possible. :math:\mathrm{eps} is too small: :math:\mathrm{eps} = \langle\mathit{\boldsymbol{value}}\rangle. The final value of :math:\mathrm{x} returned is an accurate approximation to the zero.

(errno :math:3)
The function values in the interval :math:\left[\mathrm{a}, \mathrm{b}\right] might contain a pole rather than a zero. Reducing :math:\mathrm{eps} may help in distinguishing between a pole and a zero.

.. _c05ay-py2-py-notes:

**Notes**
contfn_brent attempts to obtain an approximation to a simple zero of the function :math:f\left(x\right) given an initial interval :math:\left[a, b\right] such that :math:f\left(a\right)\times f\left(b\right)\leq 0. The same core algorithm is used by :meth:contfn_brent_rcomm whose specification should be consulted for details of the method used.

The approximation :math:x to the zero :math:\alpha is determined so that at least one of the following criteria is satisfied:

(i) :math:\left\lvert x-\alpha \right\rvert \leq \mathrm{eps},

(#) :math:\left\lvert f\left(x\right)\right\rvert \leq \mathrm{eta}.

.. _c05ay-py2-py-references:

**References**
Brent, R P, 1973, Algorithms for Minimization Without Derivatives, Prentice--Hall
"""
raise NotImplementedError

[docs]def contfn_brent_rcomm(x, y, fx, tolx, c, ind, ir=0):
r"""
contfn_brent_rcomm locates a simple zero of a continuous function in a given interval by using Brent's method, which is a combination of nonlinear interpolation, linear extrapolation and bisection.
It uses reverse communication for evaluating the function.

.. _c05az-py2-py-doc:

For full information please refer to the NAG Library document for c05az

https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05azf.html

.. _c05az-py2-py-parameters:

**Parameters**
**x** : float
On initial entry: :math:\mathrm{x} and :math:\mathrm{y} must define an initial interval :math:\left[a, b\right] containing the zero, such that :math:f\left(\mathrm{x}\right)\times f\left(\mathrm{y}\right)\leq 0.0. It is not necessary that :math:\mathrm{x} < \mathrm{y}. Otherwise, :math:\mathrm{x} and :math:\mathrm{y} must be the same entities output by a previous call to contfn_brent_rcomm.

**y** : float
On initial entry: :math:\mathrm{x} and :math:\mathrm{y} must define an initial interval :math:\left[a, b\right] containing the zero, such that :math:f\left(\mathrm{x}\right)\times f\left(\mathrm{y}\right)\leq 0.0. It is not necessary that :math:\mathrm{x} < \mathrm{y}. Otherwise, :math:\mathrm{x} and :math:\mathrm{y} must be the same entities output by a previous call to contfn_brent_rcomm.

**fx** : float
On initial entry: if :math:\mathrm{ind} = 1, :math:\mathrm{fx} need not be set.

If :math:\mathrm{ind} = -1, :math:\mathrm{fx} must contain :math:f\left(\mathrm{x}\right) for the initial value of :math:\mathrm{x}.

On intermediate entry: must contain :math:f\left(\mathrm{x}\right) for the current value of :math:\mathrm{x}.

**tolx** : float
On initial entry: the accuracy to which the zero is required. The type of error test is specified by :math:\mathrm{ir}.

**c** : float, ndarray, shape :math:\left(17\right), modified in place
On initial entry: if :math:\mathrm{ind} = 1, no elements of :math:\mathrm{c} need be set.

If :math:\mathrm{ind} = -1, :math:\mathrm{c} must contain :math:f\left(\mathrm{y}\right), other elements of :math:\mathrm{c} need not be set.

On final exit: is undefined.

**ind** : int
On initial entry: must be set to :math:1 or :math:-1.

:math:\mathrm{ind} = 1

:math:\mathrm{fx} and :math:\mathrm{c} need not be set.

:math:\mathrm{ind} = -1

:math:\mathrm{fx} and :math:\mathrm{c} must contain :math:f\left(\mathrm{x}\right) and :math:f\left(\mathrm{y}\right) respectively.

**ir** : int, optional
On initial entry: indicates the type of error test.

:math:\mathrm{ir} = 0

The test is: :math:\left\lvert \mathrm{x}-\mathrm{y}\right\rvert \leq 2.0\times \mathrm{tolx}\times \mathrm{max}\left(1.0, \left\lvert \mathrm{x}\right\rvert \right).

:math:\mathrm{ir} = 1

The test is: :math:\left\lvert \mathrm{x}-\mathrm{y}\right\rvert \leq 2.0\times \mathrm{tolx}.

:math:\mathrm{ir} = 2

The test is: :math:\left\lvert \mathrm{x}-\mathrm{y}\right\rvert \leq 2.0\times \mathrm{tolx}\times \left\lvert \mathrm{x}\right\rvert.

**Returns**
**x** : float
On intermediate exit: :math:\mathrm{x} contains the point at which :math:f must be evaluated before re-entry to the function. :math:\mathrm{y} has no external significance but must be communicated back to the next call of contfn_brent_rcomm as argument :math:\mathrm{y}.

On final exit: :math:\mathrm{x} and :math:\mathrm{y} define a smaller interval containing the zero, such that :math:f\left(\mathrm{x}\right)\times f\left(\mathrm{y}\right)\leq 0.0, and :math:\left\lvert \mathrm{x}-\mathrm{y}\right\rvert satisfies the accuracy specified by :math:\mathrm{tolx} and :math:\mathrm{ir}, unless an error has occurred. If :math:\mathrm{errno} = 4, :math:\mathrm{x} and :math:\mathrm{y} generally contain very good approximations to a pole; if :math:\mathrm{errno} = 5, :math:\mathrm{x} and :math:\mathrm{y} generally contain very good approximations to the zero (see :ref:Exceptions <c05az-py2-py-errors>). If a point :math:\mathrm{x} is found such that :math:f\left(\mathrm{x}\right) = 0.0, on final exit :math:\mathrm{x} = \mathrm{y} (in this case there is no guarantee that :math:\mathrm{x} is a simple zero). In all cases, the value returned in :math:\mathrm{x} is the better approximation to the zero.

**y** : float
On intermediate exit: :math:\mathrm{x} contains the point at which :math:f must be evaluated before re-entry to the function. :math:\mathrm{y} has no external significance but must be communicated back to the next call of contfn_brent_rcomm as argument :math:\mathrm{y}.

On final exit: :math:\mathrm{x} and :math:\mathrm{y} define a smaller interval containing the zero, such that :math:f\left(\mathrm{x}\right)\times f\left(\mathrm{y}\right)\leq 0.0, and :math:\left\lvert \mathrm{x}-\mathrm{y}\right\rvert satisfies the accuracy specified by :math:\mathrm{tolx} and :math:\mathrm{ir}, unless an error has occurred. If :math:\mathrm{errno} = 4, :math:\mathrm{x} and :math:\mathrm{y} generally contain very good approximations to a pole; if :math:\mathrm{errno} = 5, :math:\mathrm{x} and :math:\mathrm{y} generally contain very good approximations to the zero (see :ref:Exceptions <c05az-py2-py-errors>). If a point :math:\mathrm{x} is found such that :math:f\left(\mathrm{x}\right) = 0.0, on final exit :math:\mathrm{x} = \mathrm{y} (in this case there is no guarantee that :math:\mathrm{x} is a simple zero). In all cases, the value returned in :math:\mathrm{x} is the better approximation to the zero.

**ind** : int
On intermediate exit: contains :math:2, :math:3 or :math:4. The calling program must evaluate :math:f at :math:\mathrm{x}, storing the result in :math:\mathrm{fx}, and re-enter contfn_brent_rcomm with all other arguments unchanged.

On final exit: contains :math:0.

.. _c05az-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:1)
On entry, :math:f\left(\mathrm{x}\right) and :math:f\left(\mathrm{y}\right) have the same sign with neither equalling :math:0.0: :math:f\left(\mathrm{x}\right) = \langle\mathit{\boldsymbol{value}}\rangle and :math:f\left(\mathrm{y}\right) = \langle\mathit{\boldsymbol{value}}\rangle.

(errno :math:2)
On entry, :math:\mathrm{ind} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{ind} = -1, :math:1, :math:2, :math:3 or :math:4.

(errno :math:3)
On entry, :math:\mathrm{tolx} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{tolx} > 0.0.

(errno :math:3)
On entry, :math:\mathrm{ir} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{ir} = 0, :math:1 or :math:2.

**Warns**
**NagAlgorithmicWarning**
(errno :math:4)
The final interval may contain a pole rather than a zero. Note that this error exit is not completely reliable: it may be taken in extreme cases when :math:\left[\mathrm{x}, \mathrm{y}\right] contains a zero, or it may not be taken when :math:\left[\mathrm{x}, \mathrm{y}\right] contains a pole. Both these cases occur most frequently when :math:\mathrm{tolx} is large.

(errno :math:5)
The tolerance :math:\mathrm{tolx} has been set too small for the problem being solved. However, the values :math:\mathrm{x} and :math:\mathrm{y} returned may well be good approximations to the zero. :math:\mathrm{tolx} = \langle\mathit{\boldsymbol{value}}\rangle.

.. _c05az-py2-py-notes:

**Notes**
You must supply :math:\mathrm{x} and :math:\mathrm{y} to define an initial interval :math:\left[a, b\right] containing a simple zero of the function :math:f\left(x\right) (the choice of :math:\mathrm{x} and :math:\mathrm{y} must be such that :math:f\left(\mathrm{x}\right)\times f\left(\mathrm{y}\right)\leq 0.0).
The function combines the methods of bisection, nonlinear interpolation and linear extrapolation (see Dahlquist and Björck (1974)), to find a sequence of sub-intervals of the initial interval such that the final interval :math:\left[\mathrm{x}, \mathrm{y}\right] contains the zero and :math:\left\lvert \mathrm{x}-\mathrm{y}\right\rvert is less than some tolerance specified by :math:\mathrm{tolx} and :math:\mathrm{ir} (see :ref:Parameters <c05az-py2-py-parameters>).
In fact, since the intermediate intervals :math:\left[\mathrm{x}, \mathrm{y}\right] are determined only so that :math:f\left(\mathrm{x}\right)\times f\left(\mathrm{y}\right)\leq 0.0, it is possible that the final interval may contain a discontinuity or a pole of :math:f (violating the requirement that :math:f be continuous). contfn_brent_rcomm checks if the sign change is likely to correspond to a pole of :math:f and gives an error return in this case.

A feature of the algorithm used by this function is that unlike some other methods it guarantees convergence within about :math:\left(\mathrm{log}_2\left[\left(b-a\right)/\delta \right]\right)^2 function evaluations, where :math:\delta is related to the argument :math:\mathrm{tolx}.
See Brent (1973) for more details.

contfn_brent_rcomm returns to the calling program for each evaluation of :math:f\left(x\right).
On each return you should set :math:\mathrm{fx} = f\left(\mathrm{x}\right) and call contfn_brent_rcomm again.

The function is a modified version of procedure 'zeroin' given by Brent (1973).

.. _c05az-py2-py-references:

**References**
Brent, R P, 1973, Algorithms for Minimization Without Derivatives, Prentice--Hall

Bus, J C P and Dekker, T J, 1975, Two efficient algorithms with guaranteed convergence for finding a zero of a function, ACM Trans. Math. Software (1), 330--345

Dahlquist, G and Björck, Å, 1974, Numerical Methods, Prentice--Hall
"""
raise NotImplementedError

[docs]def lambertw_real(x, branch, offset):
r"""
lambertw_real returns the real values of Lambert's :math:W function :math:W\left(x\right).

.. _c05ba-py2-py-doc:

For full information please refer to the NAG Library document for c05ba

https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05baf.html

.. _c05ba-py2-py-parameters:

**Parameters**
**x** : float
If :math:\mathrm{offset} = \mathbf{True}, :math:\mathrm{x} is the offset :math:{\Delta x} from :math:{-\mathrm{exp}\left(-1\right)} of the intended argument to :math:W; that is, :math:W\left(\beta \right) is computed, where :math:\beta = {-\mathrm{exp}\left(-1\right)}+{\Delta x}.

If :math:\mathrm{offset} = \mathbf{False}, :math:\mathrm{x} is the argument :math:x of the function; that is, :math:W\left(\beta \right) is computed, where :math:\beta = x.

**branch** : int
The real branch required.

:math:\mathrm{branch} = 0

The branch :math:W_0 is selected.

:math:\mathrm{branch} = -1

The branch :math:W_{-1} is selected.

**offset** : bool
Controls whether or not :math:\mathrm{x} is being specified as an offset from :math:{-\mathrm{exp}\left(-1\right)}.

**Returns**
**lambw** : float
The real values of Lambert's :math:W function :math:W\left(x\right).

.. _c05ba-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:1)
On entry, :math:\mathrm{branch} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{branch} = 0 or :math:-1.

(errno :math:1)
On entry, :math:\mathrm{offset} = \mathbf{True} and :math:\mathrm{x} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: if :math:\mathrm{offset} = \mathbf{True} then :math:\mathrm{x}\geq 0.0.

(errno :math:1)
On entry, :math:\mathrm{offset} = \mathbf{False} and :math:\mathrm{x} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: if :math:\mathrm{offset} = \mathbf{False} then :math:\mathrm{x} \geq {-\mathrm{exp}\left({-1.0}\right)}.

(errno :math:1)
On entry, :math:\mathrm{branch} = -1, :math:\mathrm{offset} = \mathbf{True} and :math:\mathrm{x} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: if :math:\mathrm{branch} = -1 and :math:\mathrm{offset} = \mathbf{True} then :math:\mathrm{x} < {\mathrm{exp}\left({-1.0}\right)}.

(errno :math:1)
On entry, :math:\mathrm{branch} = -1, :math:\mathrm{offset} = \mathbf{False} and :math:\mathrm{x} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: if :math:\mathrm{branch} = -1 and :math:\mathrm{offset} = \mathbf{False} then :math:\mathrm{x} < 0.0.

**Warns**
**NagAlgorithmicWarning**
(errno :math:2)
For the given offset :math:\mathrm{x}, :math:W is negligibly different from :math:-1: :math:\mathrm{x} = \langle\mathit{\boldsymbol{value}}\rangle.

(errno :math:2)
:math:\mathrm{x} is close to :math:{-\mathrm{exp}\left(-1\right)}. Enter :math:\mathrm{x} as an offset to :math:{-\mathrm{exp}\left(-1\right)} for greater accuracy: :math:\mathrm{x} = \langle\mathit{\boldsymbol{value}}\rangle.

.. _c05ba-py2-py-notes:

**Notes**
lambertw_real calculates an approximate value for the real branches of Lambert's :math:W function (sometimes known as the 'product log' or 'Omega' function), which is the inverse function of

.. math::

The function :math:f is many-to-one, and so, except at :math:0, :math:W is multivalued. lambertw_real restricts :math:W and its argument :math:x to be real, resulting in a function defined for :math:x\geq {-\mathrm{exp}\left(-1\right)} and which is double valued on the interval :math:\left({-\mathrm{exp}\left(-1\right)}, 0\right).
This double-valued function is split into two real-valued branches according to the sign of :math:W\left(x\right)+1.
We denote by :math:W_0 the branch satisfying :math:W_0\left(x\right)\geq -1 for all real :math:x, and by :math:W_{-1} the branch satisfying :math:W_{-1}\left(x\right)\leq -1 for all real :math:x.
You may select your branch of interest using the argument :math:\mathrm{branch}.

The precise method used to approximate :math:W is described fully in Barry et al. (1995).
For :math:x close to :math:{-\mathrm{exp}\left(-1\right)} greater accuracy comes from evaluating :math:W\left({-\mathrm{exp}\left(-1\right)}+{\Delta x}\right) rather than :math:W\left(x\right): by setting :math:\mathrm{offset} = \mathbf{True} on entry you inform lambertw_real that you are providing :math:{\Delta x}, not :math:x, in :math:\mathrm{x}.

.. _c05ba-py2-py-references:

**References**
Barry, D J, Culligan--Hensley, P J, and Barry, S J, 1995, Real values of the :math:W -function, ACM Trans. Math. Software (21(2)), 161--171
"""
raise NotImplementedError

[docs]def lambertw_complex(branch, offset, z):
r"""
lambertw_complex computes the values of Lambert's :math:W function :math:W\left(z\right).

.. _c05bb-py2-py-doc:

For full information please refer to the NAG Library document for c05bb

https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05bbf.html

.. _c05bb-py2-py-parameters:

**Parameters**
**branch** : int
The branch required.

**offset** : bool
Controls whether or not :math:\mathrm{z} is being specified as an offset from :math:{-\mathrm{exp}\left(-1\right)}.

**z** : complex
If :math:\mathrm{offset} = \mathrm{True}, :math:\mathrm{z} is the offset :math:{\Delta z} from :math:{-\mathrm{exp}\left(-1\right)} of the intended argument to :math:W; that is, :math:W\left(\beta \right) is computed, where :math:\beta = {-\mathrm{exp}\left(-1\right)}+{\Delta z}.

If :math:\mathrm{offset} = \mathrm{False}, :math:\mathrm{z} is the argument :math:z of the function; that is, :math:W\left(\beta \right) is computed, where :math:\beta = z.

**Returns**
**w** : complex
The value :math:W\left(\beta \right): see also the description of :math:\mathrm{z}.

**resid** : float
The residual :math:\left\lvert W\left(\beta \right)\mathrm{exp}\left(W\left(\beta \right)\right)-\beta \right\rvert: see also the description of :math:\mathrm{z}.

.. _c05bb-py2-py-errors:

**Warns**
**NagAlgorithmicWarning**
(errno :math:1)
For the given offset :math:\mathrm{z}, :math:W is negligibly different from :math:-1: :math:\mathrm{Re}\left(\mathrm{z}\right) = \langle\mathit{\boldsymbol{value}}\rangle and :math:\mathrm{Im}\left(\mathrm{z}\right) = \langle\mathit{\boldsymbol{value}}\rangle.

(errno :math:1)
:math:\mathrm{z} is close to :math:{-\mathrm{exp}\left(-1\right)}. Enter :math:\mathrm{z} as an offset to :math:{-\mathrm{exp}\left(-1\right)} for greater accuracy: :math:\mathrm{Re}\left(\mathrm{z}\right) = \langle\mathit{\boldsymbol{value}}\rangle and :math:\mathrm{Im}\left(\mathrm{z}\right) = \langle\mathit{\boldsymbol{value}}\rangle.

(errno :math:2)
The iterative procedure used internally did not converge in :math:\langle\mathit{\boldsymbol{value}}\rangle iterations. Check the value of :math:\mathrm{resid} for the accuracy of :math:\mathrm{w}.

.. _c05bb-py2-py-notes:

**Notes**
lambertw_complex calculates an approximate value for Lambert's :math:W function (sometimes known as the 'product log' or 'Omega' function), which is the inverse function of

.. math::

The function :math:f is many-to-one, and so, except at :math:0, :math:W is multivalued. lambertw_complex allows you to specify the branch of :math:W on which you would like the results to lie by using the argument :math:\mathrm{branch}.
Our choice of branch cuts is as in Corless et al. (1996), and the ranges of the branches of :math:W are summarised in Figure [label omitted].

[figure omitted]

For more information about the closure of each branch, which is not displayed in Figure [label omitted], see Corless et al. (1996).
The dotted lines in the Figure denote the asymptotic boundaries of the branches, at multiples of :math:\pi.

The precise method used to approximate :math:W is as described in Corless et al. (1996).
For :math:z close to :math:{-\mathrm{exp}\left(-1\right)} greater accuracy comes from evaluating :math:W\left({-\mathrm{exp}\left(-1\right)}+{\Delta z}\right) rather than :math:W\left(z\right): by setting :math:\mathrm{offset} = \mathrm{True} on entry you inform lambertw_complex that you are providing :math:{\Delta z}, not :math:z, in :math:\mathrm{z}.

.. _c05bb-py2-py-references:

**References**
Corless, R M, Gonnet, G H, Hare, D E G, Jeffrey, D J and Knuth, D E, 1996, On the Lambert :math:W function, Advances in Comp. Math. (3), 329--359
"""
raise NotImplementedError

[docs]def sys_func_aa(fcn, x, atol=sqrt(machine.precision()), rtol=sqrt(machine.precision()), m=0, cndtol=0.0, astart=0, data=None):
r"""
sys_func_aa finds a solution of a system of nonlinear equations by fixed-point iteration using Anderson acceleration.

.. _c05mb-py2-py-doc:

For full information please refer to the NAG Library document for c05mb

https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05mbf.html

.. _c05mb-py2-py-parameters:

**Parameters**
**fcn** : callable fvec = fcn(x, data=None)
:math:\mathrm{fcn} must return the values of the functions :math:f_k at a point :math:x.

**Parameters**
**x** : float, ndarray, shape :math:\left(n\right)
The components of the point :math:x at which the functions must be evaluated.

**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.

**Returns**
**fvec** : float, array-like, shape :math:\left(n\right)
The function values :math:f_k\left(x\right).

**x** : float, array-like, shape :math:\left(n\right)
An initial guess at the solution vector, :math:\hat{x}_0.

**atol** : float, optional
The absolute convergence criterion; see :math:\mathrm{rtol}.

**rtol** : float, optional
The relative convergence criterion. At each iteration :math:\left\lVert f\left(\hat{x}_i\right)\right\rVert is computed. The iteration is deemed to have converged if :math:\left\lVert f\left(\hat{x}_i\right)\right\rVert \leq \mathrm{max}\left(\mathrm{atol}, {\mathrm{rtol}\times \left\lVert f\left(\hat{x}_0\right)\right\rVert }\right).

**m** : int, optional
:math:m, the number of previous iterates to use in Anderson acceleration. If :math:m = 0, Anderson acceleration is not used.

**cndtol** : float, optional
The maximum allowable condition number for the triangular :math:QR factor generated during Anderson acceleration. At each iteration, if the condition number exceeds :math:\mathrm{cndtol}, columns are deleted until it is sufficiently small.

If :math:\mathrm{cndtol} = 0.0, no condition number tests are performed.

**astart** : int, optional
The number of iterations by which to delay the start of Anderson acceleration.

**data** : arbitrary, optional
User-communication data for callback functions.

**Returns**
**x** : float, ndarray, shape :math:\left(n\right)
The final estimate of the solution vector.

**fvec** : float, ndarray, shape :math:\left(n\right)
The function values at the final point, :math:\mathrm{x}.

.. _c05mb-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:1)
On entry, :math:n = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:n > 0.

(errno :math:2)
On entry, :math:\mathrm{atol} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{atol}\geq 0.0.

(errno :math:3)
On entry, :math:\mathrm{rtol} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{rtol}\geq 0.0.

(errno :math:4)
On entry, :math:\mathrm{m} = \langle\mathit{\boldsymbol{value}}\rangle and :math:n = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:0\leq \mathrm{m}\leq n.

(errno :math:5)
On entry, :math:\mathrm{cndtol} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{cndtol}\geq 0.0.

(errno :math:6)
On entry, :math:\mathrm{astart} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{astart} \geq 0.

(errno :math:7)
An error occurred in evaluating the :math:QR decomposition during Anderson acceleration. This may be due to slow convergence of the iteration. Try setting the value of :math:\mathrm{cndtol}. If condition number tests are already performed, try decreasing :math:\mathrm{cndtol}.

(errno :math:11)
The iteration has diverged and subsequent iterates are too large to be computed in floating-point arithmetic.

**Warns**
**NagAlgorithmicWarning**
(errno :math:8)
The iteration is not making good progress. This failure exit may indicate that the system does not have a zero, or that the solution is very close to the origin. Rerunning sys_func_aa from a different starting point may avoid the region of difficulty.

(errno :math:9)
There have been at least :math:200\times \left(n+1\right) calls to :math:\mathrm{fcn}. Consider restarting the calculation from the point held in :math:\mathrm{x}.

**NagCallbackTerminateWarning**
(errno :math:10)
Termination requested in :math:\mathrm{fcn}.

.. _c05mb-py2-py-notes:

**Notes**
The system of equations is defined as:

.. math::
f_k\left(x_1, x_2, \ldots, x_n\right) = 0\text{, }\quad k = 1,2,\ldots,n\text{.}

This homogeneous system can readily be reformulated as

.. math::
g\left(x\right) = x\text{, }\quad x \in \mathbb{R}^n\text{.}

A standard fixed-point iteration approach is to start with an approximate solution :math:\hat{x}_0 and repeatedly apply the function :math:g until possible convergence; i.e., :math:\hat{x}_{{i+1}} = g\left(\hat{x}_i\right), until :math:\left\lVert \hat{x}_{{i+1}}-\hat{x}_i\right\rVert < \text{tol}.
Anderson acceleration uses up to :math:m previous values of :math:\hat{x} to obtain an improved estimate :math:\hat{x}_{{i+1}}.
If a standard fixed-point iteration converges, Anderson acceleration usually results in convergence in far fewer iterations (therefore, using far fewer function evaluations).

Full details of Anderson acceleration are provided in Anderson (1965).
In summary, the previous :math:m iterates are combined to form a succession of least squares problems.
These are solved using a :math:QR decomposition, which is updated at each iteration.

You are free to choose any value for :math:m, provided :math:m\leq n.
A typical choice is :math:m = 4.

Anderson acceleration is particularly useful if evaluating :math:g\left(x\right) is very expensive, in which case functions such as :meth:sys_deriv_easy or :meth:sys_func_easy, which require the Jacobian or its approximation, may perform poorly.

.. _c05mb-py2-py-references:

**References**
Anderson, D G, 1965, Iterative Procedures for Nonlinear Integral Equations, J. Assoc. Comput. Mach. (12), 547--560
"""
raise NotImplementedError

[docs]def sys_func_aa_rcomm(irevcm, x, fvec, comm, atol=sqrt(machine.precision()), rtol=sqrt(machine.precision()), m=0, cndtol=0.0, astart=0):
r"""
sys_func_aa_rcomm is a comprehensive reverse communication function that finds a solution of a system of nonlinear equations by fixed-point iteration using Anderson acceleration.

.. _c05md-py2-py-doc:

For full information please refer to the NAG Library document for c05md

https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05mdf.html

.. _c05md-py2-py-parameters:

**Parameters**
**irevcm** : int
On initial entry: must have the value :math:0.

**x** : float, ndarray, shape :math:\left(n\right), modified in place
On initial entry: an initial guess at the solution vector, :math:\hat{x}_0.

On intermediate exit: contains the current point.

On final exit: the final estimate of the solution vector.

**fvec** : float, ndarray, shape :math:\left(n\right), modified in place
On initial entry: need not be set.

On intermediate entry: if :math:\mathrm{irevcm} = 1, :math:\mathrm{fvec} must not be changed.

If :math:\mathrm{irevcm} = 2, :math:\mathrm{fvec} must be set to the values of the functions computed at the current point :math:\mathrm{x}, :math:f\left(\hat{x}_i\right).

On final exit: the function values at the final point, :math:\mathrm{x}.

**comm** : dict, communication object, modified in place
Communication structure.

On initial entry: need not be set.

**atol** : float, optional
On initial entry: the absolute convergence criterion; see :math:\mathrm{rtol}.

**rtol** : float, optional
On initial entry: the relative convergence criterion. At each iteration :math:\left\lVert f\left(\hat{x}_i\right)\right\rVert is computed. The iteration is deemed to have converged if :math:\left\lVert f\left(\hat{x}_i\right)\right\rVert \leq \mathrm{max}\left(\mathrm{atol}, {\mathrm{rtol}\times \left\lVert f\left(\hat{x}_0\right)\right\rVert }\right).

**m** : int, optional
On initial entry: :math:m, the number of previous iterates to use in Anderson acceleration. If :math:m = 0, Anderson acceleration is not used.

**cndtol** : float, optional
On initial entry: the maximum allowable condition number for the triangular :math:QR factor generated during Anderson acceleration. At each iteration, if the condition number exceeds :math:\mathrm{cndtol}, columns are deleted until it is sufficiently small.

If :math:\mathrm{cndtol} = 0.0, no condition number tests are performed.

**astart** : int, optional
On initial entry: the number of iterations by which to delay the start of Anderson acceleration.

**Returns**
**irevcm** : int
On intermediate exit: specifies what action you must take before re-entering sys_func_aa_rcomm with :math:\mathrm{irevcm} set to this value. The value of :math:\mathrm{irevcm} should be interpreted as follows:

:math:\mathrm{irevcm} = 1

Indicates the start of a new iteration. No action is required by you, but :math:\mathrm{x} and :math:\mathrm{fvec} are available for printing, and a limit on the number of iterations can be applied.

:math:\mathrm{irevcm} = 2

Indicates that before re-entry to sys_func_aa_rcomm, :math:\mathrm{fvec} must contain the function values :math:f\left(\hat{x}_i\right).

On final exit: :math:\mathrm{irevcm} = 0 and the algorithm has terminated.

.. _c05md-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:1)
On initial entry, :math:\mathrm{irevcm} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{irevcm} = 0.

(errno :math:1)
On intermediate entry, :math:\mathrm{irevcm} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{irevcm} = 1 or :math:2.

(errno :math:2)
On entry, :math:n = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:n > 0.

(errno :math:3)
On entry, :math:\mathrm{atol} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{atol}\geq 0.0.

(errno :math:4)
On entry, :math:\mathrm{rtol} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{rtol}\geq 0.0.

(errno :math:5)
On entry, :math:\mathrm{m} = \langle\mathit{\boldsymbol{value}}\rangle and :math:n = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:0\leq \mathrm{m}\leq n.

(errno :math:6)
On entry, :math:\mathrm{cndtol} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{cndtol}\geq 0.0.

(errno :math:7)
On entry, :math:\mathrm{astart} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{astart} \geq 0.

(errno :math:8)
An error occurred in evaluating the :math:QR decomposition during Anderson acceleration. This may be due to slow convergence of the iteration. Try setting the value of :math:\mathrm{cndtol}. If condition number tests are already performed, try decreasing :math:\mathrm{cndtol}.

(errno :math:9)
The iteration is not making good progress, as measured by the reduction in the norm of :math:f\left(x\right) in the last :math:\langle\mathit{\boldsymbol{value}}\rangle iterations.

(errno :math:10)
The iteration has diverged and subsequent iterates are too large to be computed in floating-point arithmetic.

.. _c05md-py2-py-notes:

**Notes**
The system of equations is defined as:

.. math::
f_k\left(x_1, x_2, \ldots, x_n\right) = 0\text{, }\quad k = 1,2,\ldots,n\text{.}

This homogeneous system can readily be reformulated as

.. math::
g\left(x\right) = x\text{, }\quad x \in \mathbb{R}^n\text{.}

A standard fixed-point iteration approach is to start with an approximate solution :math:\hat{x}_0 and repeatedly apply the function :math:g until possible convergence; i.e., :math:\hat{x}_{{i+1}} = g\left(\hat{x}_i\right), until :math:\left\lVert \hat{x}_{{i+1}}-\hat{x}_i\right\rVert < \text{tol}.
Anderson acceleration uses up to :math:m previous values of :math:\hat{x} to obtain an improved estimate :math:\hat{x}_{{i+1}}.
If a standard fixed-point iteration converges, then Anderson acceleration usually results in convergence in far fewer iterations (and, therefore, using far fewer function evaluations).

Full details of Anderson acceleration are provided in Anderson (1965).
In summary, the previous :math:m iterates are combined to form a succession of least squares problems.
These are solved using a :math:QR decomposition, which is updated at each iteration.

You are free to choose any value for :math:m, provided :math:m\leq n.
A typical choice is :math:m = 4.

Anderson acceleration is particularly useful if evaluating :math:g\left(x\right) is very expensive, in which case functions such as :meth:sys_deriv_rcomm or :meth:sys_func_rcomm, which require the Jacobian or its approximation, may perform poorly.

.. _c05md-py2-py-references:

**References**
Anderson, D G, 1965, Iterative Procedures for Nonlinear Integral Equations, J. Assoc. Comput. Mach. (12), 547--560
"""
raise NotImplementedError

[docs]def sys_func_easy(fcn, x, xtol=sqrt(machine.precision()), data=None):
r"""
sys_func_easy is an easy-to-use function that finds a solution of a system of nonlinear equations by a modification of the Powell hybrid method.

.. _c05qb-py2-py-doc:

For full information please refer to the NAG Library document for c05qb

https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05qbf.html

.. _c05qb-py2-py-parameters:

**Parameters**
**fcn** : callable fvec = fcn(x, data=None)
:math:\mathrm{fcn} must return the values of the functions :math:f_i at a point :math:x.

**Parameters**
**x** : float, ndarray, shape :math:\left(n\right)
The components of the point :math:x at which the functions must be evaluated.

**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.

**Returns**
**fvec** : float, array-like, shape :math:\left(n\right)
The function values :math:f_i\left(x\right).

**x** : float, array-like, shape :math:\left(n\right)
An initial guess at the solution vector.

**xtol** : float, optional
The accuracy in :math:\mathrm{x} to which the solution is required.

**data** : arbitrary, optional
User-communication data for callback functions.

**Returns**
**x** : float, ndarray, shape :math:\left(n\right)
The final estimate of the solution vector.

**fvec** : float, ndarray, shape :math:\left(n\right)
The function values at the final point returned in :math:\mathrm{x}.

.. _c05qb-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:11)
On entry, :math:n = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:n > 0.

(errno :math:12)
On entry, :math:\mathrm{xtol} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{xtol}\geq 0.0.

**Warns**
**NagAlgorithmicWarning**
(errno :math:2)
There have been at least :math:200\times \left(n+1\right) calls to :math:\mathrm{fcn}. Consider restarting the calculation from the point held in :math:\mathrm{x}.

(errno :math:3)
No further improvement in the solution is possible. :math:\mathrm{xtol} is too small: :math:\mathrm{xtol} = \langle\mathit{\boldsymbol{value}}\rangle.

(errno :math:4)
The iteration is not making good progress. This failure exit may indicate that the system does not have a zero, or that the solution is very close to the origin (see Accuracy <https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05qbf.html#accuracy>__). Otherwise, rerunning sys_func_easy from a different starting point may avoid the region of difficulty.

**NagCallbackTerminateWarning**
(errno :math:5)
Termination requested in :math:\mathrm{fcn}.

.. _c05qb-py2-py-notes:

**Notes**
The system of equations is defined as:

.. math::
f_i\left(x_1, x_2, \ldots, x_n\right) = 0\text{, }i = 1,2,\ldots,n\text{.}

sys_func_easy is based on the MINPACK routine HYBRD1 (see Moré et al. (1980)).
It chooses the correction at each step as a convex combination of the Newton and scaled gradient directions.
The Jacobian is updated by the rank-1 method of Broyden.
At the starting point, the Jacobian is approximated by forward differences, but these are not used again until the rank-1 method fails to produce satisfactory progress.
For more details see Powell (1970).

.. _c05qb-py2-py-references:

**References**
Moré, J J, Garbow, B S and Hillstrom, K E, 1980, User guide for MINPACK-1, Technical Report ANL-80-74, Argonne National Laboratory

Powell, M J D, 1970, A hybrid method for nonlinear algebraic equations, Numerical Methods for Nonlinear Algebraic Equations, (ed P Rabinowitz), Gordon and Breach
"""
raise NotImplementedError

[docs]def sys_func_expert(fcnval, x, ml, mu, mode, diag, nprint, fcnmon=None, xtol=sqrt(machine.precision()), maxfev=None, epsfcn=0.0, factor=100.0, data=None):
r"""
sys_func_expert is a comprehensive function that finds a solution of a system of nonlinear equations by a modification of the Powell hybrid method.

.. _c05qc-py2-py-doc:

For full information please refer to the NAG Library document for c05qc

https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05qcf.html

.. _c05qc-py2-py-parameters:

**Parameters**
**fcnval** : callable fvec = fcnval(x, data=None)
:math:\mathrm{fcnval} must return the values of the functions :math:f_i at a point :math:x.

**Parameters**
**x** : float, ndarray, shape :math:\left(n\right)
The components of the point :math:x at which the functions must be evaluated.

**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.

**Returns**
**fvec** : float, array-like, shape :math:\left(n\right)
:math:\mathrm{fvec} must contain the function values :math:f_i\left(x\right).

**x** : float, array-like, shape :math:\left(n\right)
An initial guess at the solution vector.

**ml** : int
The number of subdiagonals within the band of the Jacobian matrix. (If the Jacobian is not banded, or you are unsure, set :math:\mathrm{ml} = n-1.)

**mu** : int
The number of superdiagonals within the band of the Jacobian matrix. (If the Jacobian is not banded, or you are unsure, set :math:\mathrm{mu} = n-1.)

**mode** : int
Indicates whether or not you have provided scaling factors in :math:\mathrm{diag}.

If :math:\mathrm{mode} = 2, the scaling must have been specified in :math:\mathrm{diag}.

Otherwise, if :math:\mathrm{mode} = 1, the variables will be scaled internally.

**diag** : float, array-like, shape :math:\left(n\right)
If :math:\mathrm{mode} = 2, :math:\mathrm{diag} must contain multiplicative scale factors for the variables.

If :math:\mathrm{mode} = 1, :math:\mathrm{diag} need not be set.

**nprint** : int
Indicates whether (and how often) calls to :math:\mathrm{fcnmon} are to be made.

:math:\mathrm{nprint}\leq 0

:math:\mathrm{nprint} > 0

:math:\mathrm{fcnmon} is called at the beginning of the first iteration, every :math:\mathrm{nprint} iterations thereafter and immediately before the return from sys_func_expert.

**fcnmon** : None or callable fcnmon(x, fvec, data=None), optional
Note: if this argument is **None** then a NAG-supplied facility will be used.

:math:\mathrm{fcnmon} may be used to monitor the progress of the algorithm.

**Parameters**
**x** : float, ndarray, shape :math:\left(n\right)
The components of the current evaluation point :math:x.

**fvec** : float, ndarray, shape :math:\left(n\right)
:math:\mathrm{fvec} contains the function values :math:f_i\left(x\right).

**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.

**xtol** : float, optional
The accuracy in :math:\mathrm{x} to which the solution is required.

**maxfev** : None or int, optional
Note: if this argument is **None** then a default value will be used, determined as follows: :math:200\times \left(n+1\right).

The maximum number of calls to :math:\mathrm{fcnval}. sys_func_expert will exit with :math:\mathrm{errno} = 2, if, at the end of an iteration, the number of calls to :math:\mathrm{fcnval} exceeds :math:\mathrm{maxfev}.

**epsfcn** : float, optional
A rough estimate of the largest relative error in the functions. It is used in determining a suitable step for a forward difference approximation to the Jacobian. If :math:\mathrm{epsfcn} is less than machine precision (returned by :meth:machine.precision <naginterfaces.library.machine.precision>) then machine precision is used. Consequently a value of :math:0.0 will often be suitable.

**factor** : float, optional
A quantity to be used in determining the initial step bound. In most cases, :math:\mathrm{factor} should lie between :math:0.1 and :math:100.0. (The step bound is :math:\mathrm{factor}\times \left\lVert \mathrm{diag}\times \mathrm{x}\right\rVert_2 if this is nonzero; otherwise the bound is :math:\mathrm{factor}.)

**data** : arbitrary, optional
User-communication data for callback functions.

**Returns**
**x** : float, ndarray, shape :math:\left(n\right)
The final estimate of the solution vector.

**fvec** : float, ndarray, shape :math:\left(n\right)
The function values at the final point returned in :math:\mathrm{x}.

**diag** : float, ndarray, shape :math:\left(n\right)
The scale factors actually used (computed internally if :math:\mathrm{mode} = 1).

**nfev** : int
The number of calls made to :math:\mathrm{fcnval}.

**fjac** : float, ndarray, shape :math:\left(n, n\right)
The orthogonal matrix :math:Q produced by the :math:QR factorization of the final approximate Jacobian.

**r** : float, ndarray, shape :math:\left(n\times \left(n+1\right)/2\right)
The upper triangular matrix :math:R produced by the :math:QR factorization of the final approximate Jacobian, stored row-wise.

**qtf** : float, ndarray, shape :math:\left(n\right)
The vector :math:Q^\mathrm{T}f.

.. _c05qc-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:11)
On entry, :math:n = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:n > 0.

(errno :math:12)
On entry, :math:\mathrm{xtol} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{xtol}\geq 0.0.

(errno :math:13)
On entry, :math:\mathrm{mode} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{mode} = 1 or :math:2.

(errno :math:14)
On entry, :math:\mathrm{factor} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{factor} > 0.0.

(errno :math:15)
On entry, :math:\mathrm{mode} = 2 and :math:\mathrm{diag} contained a non-positive element.

(errno :math:16)
On entry, :math:\mathrm{ml} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{ml} \geq 0.

(errno :math:17)
On entry, :math:\mathrm{mu} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{mu} \geq 0.

(errno :math:18)
On entry, :math:\mathrm{maxfev} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{maxfev} > 0.

**Warns**
**NagAlgorithmicWarning**
(errno :math:2)
There have been at least :math:\mathrm{maxfev} calls to :math:\mathrm{fcnval}: :math:\mathrm{maxfev} = \langle\mathit{\boldsymbol{value}}\rangle. Consider restarting the calculation from the final point held in :math:\mathrm{x}.

(errno :math:3)
No further improvement in the solution is possible. :math:\mathrm{xtol} is too small: :math:\mathrm{xtol} = \langle\mathit{\boldsymbol{value}}\rangle.

(errno :math:4)
The iteration is not making good progress, as measured by the improvement from the last :math:\langle\mathit{\boldsymbol{value}}\rangle Jacobian evaluations.

(errno :math:5)
The iteration is not making good progress, as measured by the improvement from the last :math:\langle\mathit{\boldsymbol{value}}\rangle iterations.

**NagCallbackTerminateWarning**
(errno :math:6)
Termination requested in :math:\mathrm{fcnval} or :math:\mathrm{fcnmon}.

.. _c05qc-py2-py-notes:

**Notes**
The system of equations is defined as:

.. math::
f_i\left(x_1, x_2, \ldots, x_n\right) = 0\text{, }i = 1,2,\ldots,n\text{.}

sys_func_expert is based on the MINPACK routine HYBRD (see Moré et al. (1980)).
It chooses the correction at each step as a convex combination of the Newton and scaled gradient directions.
The Jacobian is updated by the rank-1 method of Broyden.
At the starting point, the Jacobian is approximated by forward differences, but these are not used again until the rank-1 method fails to produce satisfactory progress.
For more details see Powell (1970).

.. _c05qc-py2-py-references:

**References**
Moré, J J, Garbow, B S and Hillstrom, K E, 1980, User guide for MINPACK-1, Technical Report ANL-80-74, Argonne National Laboratory

Powell, M J D, 1970, A hybrid method for nonlinear algebraic equations, Numerical Methods for Nonlinear Algebraic Equations, (ed P Rabinowitz), Gordon and Breach
"""
raise NotImplementedError

[docs]def sys_func_rcomm(irevcm, x, fvec, ml, mu, mode, diag, fjac, r, qtf, comm, xtol=sqrt(machine.precision()), epsfcn=0.0, factor=100.0):
r"""
sys_func_rcomm is a comprehensive reverse communication function that finds a solution of a system of nonlinear equations by a modification of the Powell hybrid method.

.. _c05qd-py2-py-doc:

For full information please refer to the NAG Library document for c05qd

https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05qdf.html

.. _c05qd-py2-py-parameters:

**Parameters**
**irevcm** : int
On initial entry: must have the value :math:0.

**x** : float, ndarray, shape :math:\left(n\right), modified in place
On initial entry: an initial guess at the solution vector.

On intermediate exit: contains the current point.

On final exit: the final estimate of the solution vector.

**fvec** : float, ndarray, shape :math:\left(n\right), modified in place
On initial entry: need not be set.

On intermediate entry: if :math:\mathrm{irevcm} = 1, :math:\mathrm{fvec} must not be changed.

If :math:\mathrm{irevcm} = 2, :math:\mathrm{fvec} must be set to the values of the functions computed at the current point :math:\mathrm{x}.

On final exit: the function values at the final point, :math:\mathrm{x}.

**ml** : int
On initial entry: the number of subdiagonals within the band of the Jacobian matrix. (If the Jacobian is not banded, or you are unsure, set :math:\mathrm{ml} = n-1.)

**mu** : int
On initial entry: the number of superdiagonals within the band of the Jacobian matrix. (If the Jacobian is not banded, or you are unsure, set :math:\mathrm{mu} = n-1.)

**mode** : int
On initial entry: indicates whether or not you have provided scaling factors in :math:\mathrm{diag}.

If :math:\mathrm{mode} = 2, the scaling must have been supplied in :math:\mathrm{diag}.

Otherwise, if :math:\mathrm{mode} = 1, the variables will be scaled internally.

**diag** : float, ndarray, shape :math:\left(n\right), modified in place
On entry: if :math:\mathrm{mode} = 2, :math:\mathrm{diag} must contain multiplicative scale factors for the variables.

If :math:\mathrm{mode} = 1, :math:\mathrm{diag} need not be set.

On exit: the scale factors actually used (computed internally if :math:\mathrm{mode} = 1).

**fjac** : float, ndarray, shape :math:\left(n, n\right), modified in place
On initial entry: need not be set.

On intermediate exit: must not be changed.

On final exit: the orthogonal matrix :math:Q produced by the :math:QR factorization of the final approximate Jacobian.

**r** : float, ndarray, shape :math:\left(n\times \left(n+1\right)/2\right), modified in place
On initial entry: need not be set.

On intermediate exit: must not be changed.

On final exit: the upper triangular matrix :math:R produced by the :math:QR factorization of the final approximate Jacobian, stored row-wise.

**qtf** : float, ndarray, shape :math:\left(n\right), modified in place
On initial entry: need not be set.

On intermediate exit: must not be changed.

On final exit: the vector :math:Q^\mathrm{T}f.

**comm** : dict, communication object, modified in place
Communication structure.

On initial entry: need not be set.

**xtol** : float, optional
On initial entry: the accuracy in :math:\mathrm{x} to which the solution is required.

**epsfcn** : float, optional
On initial entry: the order of the largest relative error in the functions. It is used in determining a suitable step for a forward difference approximation to the Jacobian. If :math:\mathrm{epsfcn} is less than machine precision (returned by :meth:machine.precision <naginterfaces.library.machine.precision>) then machine precision is used. Consequently a value of :math:0.0 will often be suitable.

**factor** : float, optional
On initial entry: a quantity to be used in determining the initial step bound. In most cases, :math:\mathrm{factor} should lie between :math:0.1 and :math:100.0. (The step bound is :math:\mathrm{factor}\times \left\lVert \mathrm{diag}\times \mathrm{x}\right\rVert_2 if this is nonzero; otherwise the bound is :math:\mathrm{factor}.)

**Returns**
**irevcm** : int
On intermediate exit: specifies what action you must take before re-entering sys_func_rcomm with :math:\mathrm{irevcm} set to this value. The value of :math:\mathrm{irevcm} should be interpreted as follows:

:math:\mathrm{irevcm} = 1

Indicates the start of a new iteration. No action is required by you, but :math:\mathrm{x} and :math:\mathrm{fvec} are available for printing.

:math:\mathrm{irevcm} = 2

Indicates that before re-entry to sys_func_rcomm, :math:\mathrm{fvec} must contain the function values :math:f_i\left(x\right).

On final exit: :math:\mathrm{irevcm} = 0 and the algorithm has terminated.

.. _c05qd-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:2)
On entry, :math:\mathrm{irevcm} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{irevcm} = 0, :math:1 or :math:2.

(errno :math:11)
On entry, :math:n = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:n > 0.

(errno :math:12)
On entry, :math:\mathrm{xtol} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{xtol}\geq 0.0.

(errno :math:13)
On entry, :math:\mathrm{mode} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{mode} = 1 or :math:2.

(errno :math:14)
On entry, :math:\mathrm{factor} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{factor} > 0.0.

(errno :math:15)
On entry, :math:\mathrm{mode} = 2 and :math:\mathrm{diag} contained a non-positive element.

(errno :math:16)
On entry, :math:\mathrm{ml} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{ml} \geq 0.

(errno :math:17)
On entry, :math:\mathrm{mu} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{mu} \geq 0.

**Warns**
**NagAlgorithmicWarning**
(errno :math:3)
No further improvement in the solution is possible. :math:\mathrm{xtol} is too small: :math:\mathrm{xtol} = \langle\mathit{\boldsymbol{value}}\rangle.

(errno :math:4)
The iteration is not making good progress, as measured by the improvement from the last :math:\langle\mathit{\boldsymbol{value}}\rangle Jacobian evaluations.

(errno :math:5)
The iteration is not making good progress, as measured by the improvement from the last :math:\langle\mathit{\boldsymbol{value}}\rangle iterations.

.. _c05qd-py2-py-notes:

**Notes**
The system of equations is defined as:

.. math::
f_i\left(x_1, x_2, \ldots, x_n\right) = 0\text{, }\quad i = 1,2,\ldots,n\text{.}

sys_func_rcomm is based on the MINPACK routine HYBRD (see Moré et al. (1980)).
It chooses the correction at each step as a convex combination of the Newton and scaled gradient directions.
The Jacobian is updated by the rank-1 method of Broyden.
At the starting point, the Jacobian is approximated by forward differences, but these are not used again until the rank-1 method fails to produce satisfactory progress.
For more details see Powell (1970).

.. _c05qd-py2-py-references:

**References**
Moré, J J, Garbow, B S and Hillstrom, K E, 1980, User guide for MINPACK-1, Technical Report ANL-80-74, Argonne National Laboratory

Powell, M J D, 1970, A hybrid method for nonlinear algebraic equations, Numerical Methods for Nonlinear Algebraic Equations, (ed P Rabinowitz), Gordon and Breach

--------
:meth:naginterfaces.library.examples.roots.sys_func_rcomm_ex.main
"""
raise NotImplementedError

[docs]def sparsys_func_easy(fcn, x, xtol=sqrt(machine.precision()), nnz=None, comm=None, data=None):
r"""
sparsys_func_easy is an easy-to-use function that finds a solution of a sparse system of nonlinear equations by a modification of the Powell hybrid method.

.. _c05qs-py2-py-doc:

For full information please refer to the NAG Library document for c05qs

https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05qsf.html

.. _c05qs-py2-py-parameters:

**Parameters**
**fcn** : callable fcn(indf, x, fvec, data=None)
:math:\mathrm{fcn} must return the values of the functions :math:f_i at a point :math:x.

**Parameters**
**indf** : int, ndarray, shape :math:\left(\textit{lindf}\right)
Note: this argument represents array indices; the values supplied will be base-1.

:math:\mathrm{indf} specifies the indices :math:i for which values of :math:f_i\left(x\right) must be computed. The indices are specified in strictly ascending order.

**x** : float, ndarray, shape :math:\left(n\right)
The components of the point :math:x at which the functions must be evaluated. :math:\mathrm{x}[i-1] contains the coordinate :math:x_i.

**fvec** : float, ndarray, shape :math:\left(n\right), to be modified in place
On exit: :math:\mathrm{fvec}[i-1] must contain the function values :math:f_i\left(x\right), for all indices :math:i in :math:\mathrm{indf}.

**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.

**x** : float, array-like, shape :math:\left(n\right)
An initial guess at the solution vector. :math:\mathrm{x}[i-1] must contain the coordinate :math:x_i.

**xtol** : float, optional
The accuracy in :math:\mathrm{x} to which the solution is required.

**nnz** : None or int, optional
Note: if this argument is **None** then a default value will be used, determined as follows: :math:n\times n.

If supplied :math:\mathrm{nnz} should be an upper estimate of the number of nonzero entries found in the Jacobian.

The size of the data that is self communicated to this routine via :math:\mathrm{comm} is a function of :math:\mathrm{nnz}.

A suitable overestimate for :math:\mathrm{nnz} is :math:n^2, but this may be suboptimal when :math:\textit{n} is large.

If an existing sparsity pattern for the problem is communicated via :math:\mathrm{comm}, :math:\mathrm{nnz} will be ignored.

The value that was supplied in the call to sparsys_func_easy which computed that sparsity pattern will be extracted from :math:\mathrm{comm} and used instead.

**comm** : None or dict, communication object, optional, modified in place
Optionally, on entry: an existing sparsity pattern for the Jacobian, computed by a subsequent call. Otherwise, supplying an uninitialized object indicates it is the first time sparsys_func_easy is called for this specific problem.

On exit, if not **None** on entry: will be populated with the sparsity pattern of the Jacobian. On subsequent calls, the same communication object can be passed in again if the problem has a Jacobian of the same sparsity pattern.

**data** : arbitrary, optional
User-communication data for callback functions.

**Returns**
**x** : float, ndarray, shape :math:\left(n\right)
The final estimate of the solution vector.

**fvec** : float, ndarray, shape :math:\left(n\right)
The function values at the final point returned in :math:\mathrm{x}. :math:\mathrm{fvec}[i-1] contains the function values :math:f_i.

**actual_nnz** : int
The number of nonzero entries found in the Jacobian.

.. _c05qs-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:6)
On entry, the supplied :math:\mathrm{nnz} = \langle\mathit{\boldsymbol{value}}\rangle was too small for this problem.

(errno :math:7)
On entry, the supplied :math:\mathrm{nnz} = \langle\mathit{\boldsymbol{value}}\rangle was too small for this problem.

(errno :math:11)
On entry, :math:n = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:n > 0.

(errno :math:12)
On entry, :math:\mathrm{xtol} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{xtol}\geq 0.0.

**Warns**
**NagAlgorithmicWarning**
(errno :math:2)
There have been at least :math:200\times \left(n+1\right) calls to :math:\mathrm{fcn}. Consider restarting the calculation from the point held in :math:\mathrm{x} and supplying the populated sparsity pattern from :math:\mathrm{comm}.

(errno :math:3)
No further improvement in the solution is possible. :math:\mathrm{xtol} is too small: :math:\mathrm{xtol} = \langle\mathit{\boldsymbol{value}}\rangle.

(errno :math:4)
The iteration is not making good progress. This failure exit may indicate that the system does not have a zero, or that the solution is very close to the origin (see Accuracy <https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05qsf.html#accuracy>__). Otherwise, rerunning sparsys_func_easy from a different starting point may avoid the region of difficulty. The condition number of the Jacobian is :math:\langle\mathit{\boldsymbol{value}}\rangle.

**NagCallbackTerminateWarning**
(errno :math:5)
Termination requested in :math:\mathrm{fcn}.

.. _c05qs-py2-py-notes:

**Notes**
The system of equations is defined as:

.. math::
f_i\left(x_1, x_2, \ldots, x_n\right) = 0\text{, }i = 1,2,\ldots,n\text{.}

sparsys_func_easy is based on the MINPACK routine HYBRD1 (see Moré et al. (1980)).
It chooses the correction at each step as a convex combination of the Newton and scaled gradient directions.
The Jacobian is updated by the sparse rank-1 method of Schubert (see Schubert (1970)).
At the starting point, the sparsity pattern is determined and the Jacobian is approximated by forward differences, but these are not used again until the rank-1 method fails to produce satisfactory progress.
Then, the sparsity structure is used to recompute an approximation to the Jacobian by forward differences with the least number of function evaluations.
The function you supply must be able to compute only the requested subset of the function values.
The sparse Jacobian linear system is solved at each iteration with :meth:sparse.direct_real_gen_lu <naginterfaces.library.sparse.direct_real_gen_lu> computing the Newton step.
For more details see Powell (1970) and Broyden (1965).

.. _c05qs-py2-py-references:

**References**
Broyden, C G, 1965, A class of methods for solving nonlinear simultaneous equations, Mathematics of Computation (19(92)), 577--593

Moré, J J, Garbow, B S and Hillstrom, K E, 1980, User guide for MINPACK-1, Technical Report ANL-80-74, Argonne National Laboratory

Powell, M J D, 1970, A hybrid method for nonlinear algebraic equations, Numerical Methods for Nonlinear Algebraic Equations, (ed P Rabinowitz), Gordon and Breach

Schubert, L K, 1970, Modification of a quasi-Newton method for nonlinear equations with a sparse Jacobian, Mathematics of Computation (24(109)), 27--30
"""
raise NotImplementedError

[docs]def sys_deriv_easy(fcnval, fcnjac, x, xtol=sqrt(machine.precision()), data=None, spiked_sorder='C'):
r"""
sys_deriv_easy is an easy-to-use function that finds a solution of a system of nonlinear equations by a modification of the Powell hybrid method.
You must provide the Jacobian.

.. _c05rb-py2-py-doc:

For full information please refer to the NAG Library document for c05rb

https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05rbf.html

.. _c05rb-py2-py-parameters:

**Parameters**
**fcnval** : callable fvec = fcnval(x, data=None)
:math:\mathrm{fcnval} must return the values of the functions :math:f_i at a point :math:x.

**Parameters**
**x** : float, ndarray, shape :math:\left(n\right)
The components of the point :math:x at which the functions must be evaluated.

**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.

**Returns**
**fvec** : float, array-like, shape :math:\left(n\right)
:math:\mathrm{fvec} must contain the function values :math:f_i\left(x\right).

**fcnjac** : callable fjac = fcnjac(x, data=None)
:math:\mathrm{fcnjac} must return the Jacobian at :math:x.

**Parameters**
**x** : float, ndarray, shape :math:\left(n\right)
The components of the point :math:x at which the Jacobian must be evaluated.

**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.

**Returns**
**fjac** : float, array-like, shape :math:\left(n, n\right)
:math:\mathrm{fjac}[\textit{i}-1,\textit{j}-1] must contain the value of :math:\frac{{\partial f_{\textit{i}}}}{{\partial x_{\textit{j}}}} at the point :math:x, for :math:\textit{j} = 1,2,\ldots,n, for :math:\textit{i} = 1,2,\ldots,n.

**x** : float, array-like, shape :math:\left(n\right)
An initial guess at the solution vector.

**xtol** : float, optional
The accuracy in :math:\mathrm{x} to which the solution is required.

**data** : arbitrary, optional
User-communication data for callback functions.

**spiked_sorder** : str, optional
If :math:\mathrm{fjac} in :math:\mathrm{fcnjac} is spiked (i.e., has unit extent in all but one dimension, or has size :math:1), :math:\mathrm{spiked\_sorder} selects the storage order to associate with it in the NAG Engine:

spiked_sorder = :math:\texttt{'C'}
row-major storage will be used;

spiked_sorder = :math:\texttt{'F'}
column-major storage will be used.

**Returns**
**x** : float, ndarray, shape :math:\left(n\right)
The final estimate of the solution vector.

**fvec** : float, ndarray, shape :math:\left(n\right)
The function values at the final point returned in :math:\mathrm{x}.

**fjac** : float, ndarray, shape :math:\left(n, n\right)
The orthogonal matrix :math:Q produced by the :math:QR factorization of the final approximate Jacobian.

.. _c05rb-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:11)
On entry, :math:n = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:n > 0.

(errno :math:12)
On entry, :math:\mathrm{xtol} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{xtol}\geq 0.0.

**Warns**
**NagAlgorithmicWarning**
(errno :math:2)
There have been at least :math:100\times \left(n+1\right) calls to :math:\mathrm{fcnval} and :math:\mathrm{fcnjac} combined. Consider restarting the calculation from the point held in :math:\mathrm{x}.

(errno :math:3)
No further improvement in the solution is possible. :math:\mathrm{xtol} is too small: :math:\mathrm{xtol} = \langle\mathit{\boldsymbol{value}}\rangle.

(errno :math:4)
The iteration is not making good progress. This failure exit may indicate that the system does not have a zero, or that the solution is very close to the origin (see Accuracy <https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05rbf.html#accuracy>__). Otherwise, rerunning sys_deriv_easy from a different starting point may avoid the region of difficulty.

**NagCallbackTerminateWarning**
(errno :math:5)
Termination requested in :math:\mathrm{fcnval} or :math:\mathrm{fcnjac}.

.. _c05rb-py2-py-notes:

**Notes**
The system of equations is defined as:

.. math::
f_i\left(x_1, x_2, \ldots, x_n\right) = 0\text{, }\quad i = 1,2,\ldots,n\text{.}

sys_deriv_easy is based on the MINPACK routine HYBRJ1 (see Moré et al. (1980)).
It chooses the correction at each step as a convex combination of the Newton and scaled gradient directions.
The Jacobian is updated by the rank-1 method of Broyden.
At the starting point, the Jacobian is requested, but it is not asked for again until the rank-1 method fails to produce satisfactory progress.
For more details see Powell (1970).

.. _c05rb-py2-py-references:

**References**
Moré, J J, Garbow, B S and Hillstrom, K E, 1980, User guide for MINPACK-1, Technical Report ANL-80-74, Argonne National Laboratory

Powell, M J D, 1970, A hybrid method for nonlinear algebraic equations, Numerical Methods for Nonlinear Algebraic Equations, (ed P Rabinowitz), Gordon and Breach
"""
raise NotImplementedError

[docs]def sys_deriv_expert(fcnval, fcnjac, fcnmon, x, mode, diag, nprint, xtol=sqrt(machine.precision()), maxfev=None, factor=100.0, data=None, spiked_sorder='C'):
r"""
sys_deriv_expert is a comprehensive function that finds a solution of a system of nonlinear equations by a modification of the Powell hybrid method.
You must provide the Jacobian.

.. _c05rc-py2-py-doc:

For full information please refer to the NAG Library document for c05rc

https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05rcf.html

.. _c05rc-py2-py-parameters:

**Parameters**
**fcnval** : callable fvec = fcnval(x, data=None)
:math:\mathrm{fcnval} must return the values of the functions :math:f_i at a point :math:x.

**Parameters**
**x** : float, ndarray, shape :math:\left(n\right)
The components of the point :math:x at which the functions must be evaluated.

**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.

**Returns**
**fvec** : float, array-like, shape :math:\left(n\right)
:math:\mathrm{fvec} must contain the function values :math:f_i\left(x\right).

**fcnjac** : callable fjac = fcnjac(x, data=None)
:math:\mathrm{fcnjac} must return the Jacobian at :math:x.

**Parameters**
**x** : float, ndarray, shape :math:\left(n\right)
The components of the point :math:x at which the Jacobian must be evaluated.

**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.

**Returns**
**fjac** : float, array-like, shape :math:\left(n, n\right)
:math:\mathrm{fjac}[\textit{i}-1,\textit{j}-1] must contain the value of :math:\frac{{\partial f_{\textit{i}}}}{{\partial x_{\textit{j}}}} at the point :math:x, for :math:\textit{j} = 1,2,\ldots,n, for :math:\textit{i} = 1,2,\ldots,n.

**fcnmon** : callable fcnmon(x, fvec, fjac, data=None)
:math:\mathrm{fcnmon} may be used to monitor the progress of the algorithm.

**Parameters**
**x** : float, ndarray, shape :math:\left(n\right)
The components of the current evaluation point :math:x.

**fvec** : float, ndarray, shape :math:\left(n\right)
:math:\mathrm{fvec} contains the function values :math:f_i\left(x\right).

**fjac** : float, ndarray, shape :math:\left(n, n\right)
:math:\mathrm{fjac}[\textit{i}-1,\textit{j}-1] contains the value of :math:\frac{{\partial f_{\textit{i}}}}{{\partial x_{\textit{j}}}} at the point :math:x, for :math:\textit{j} = 1,2,\ldots,n, for :math:\textit{i} = 1,2,\ldots,n.

**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.

**x** : float, array-like, shape :math:\left(n\right)
An initial guess at the solution vector.

**mode** : int
Indicates whether or not you have provided scaling factors in :math:\mathrm{diag}.

If :math:\mathrm{mode} = 2, the scaling must have been specified in :math:\mathrm{diag}.

Otherwise, if :math:\mathrm{mode} = 1, the variables will be scaled internally.

**diag** : float, array-like, shape :math:\left(n\right)
If :math:\mathrm{mode} = 2, :math:\mathrm{diag} must contain multiplicative scale factors for the variables.

If :math:\mathrm{mode} = 1, :math:\mathrm{diag} need not be set.

**nprint** : int
Indicates whether (and how often) calls to :math:\mathrm{fcnmon} are to be made.

:math:\mathrm{nprint}\leq 0

:math:\mathrm{nprint} > 0

:math:\mathrm{fcnmon} is called at the beginning of the first iteration, every :math:\mathrm{nprint} iterations thereafter and immediately before the return from sys_deriv_expert.

**xtol** : float, optional
The accuracy in :math:\mathrm{x} to which the solution is required.

**maxfev** : None or int, optional
Note: if this argument is **None** then a default value will be used, determined as follows: :math:100\times \left(n+1\right).

The maximum number of calls to :math:\mathrm{fcnval} and :math:\mathrm{fcnjac} combined. sys_deriv_expert will exit with :math:\mathrm{errno} = 2, if, at the end of an iteration, the number of calls to :math:\mathrm{fcnval} and :math:\mathrm{fcnjac} combined exceeds :math:\mathrm{maxfev}.

**factor** : float, optional
A quantity to be used in determining the initial step bound. In most cases, :math:\mathrm{factor} should lie between :math:0.1 and :math:100.0. (The step bound is :math:\mathrm{factor}\times \left\lVert \mathrm{diag}\times \mathrm{x}\right\rVert_2 if this is nonzero; otherwise the bound is :math:\mathrm{factor}.)

**data** : arbitrary, optional
User-communication data for callback functions.

**spiked_sorder** : str, optional
If :math:\mathrm{fjac} in :math:\mathrm{fcnjac} is spiked (i.e., has unit extent in all but one dimension, or has size :math:1), :math:\mathrm{spiked\_sorder} selects the storage order to associate with it in the NAG Engine:

spiked_sorder = :math:\texttt{'C'}
row-major storage will be used;

spiked_sorder = :math:\texttt{'F'}
column-major storage will be used.

**Returns**
**x** : float, ndarray, shape :math:\left(n\right)
The final estimate of the solution vector.

**fvec** : float, ndarray, shape :math:\left(n\right)
The function values at the final point returned in :math:\mathrm{x}.

**fjac** : float, ndarray, shape :math:\left(n, n\right)
The orthogonal matrix :math:Q produced by the :math:QR factorization of the final approximate Jacobian.

**diag** : float, ndarray, shape :math:\left(n\right)
The scale factors actually used (computed internally if :math:\mathrm{mode} = 1).

**nfev** : int
The number of calls made to :math:\mathrm{fcnval} and :math:\mathrm{fcnjac} combined to evaluate the functions.

**njev** : int
The number of calls made to :math:\mathrm{fcnjac}.

**r** : float, ndarray, shape :math:\left(n\times \left(n+1\right)/2\right)
The upper triangular matrix :math:R produced by the :math:QR factorization of the final approximate Jacobian, stored row-wise.

**qtf** : float, ndarray, shape :math:\left(n\right)
The vector :math:Q^\mathrm{T}f.

.. _c05rc-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:11)
On entry, :math:n = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:n > 0.

(errno :math:12)
On entry, :math:\mathrm{xtol} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{xtol}\geq 0.0.

(errno :math:13)
On entry, :math:\mathrm{mode} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{mode} = 1 or :math:2.

(errno :math:14)
On entry, :math:\mathrm{factor} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{factor} > 0.0.

(errno :math:15)
On entry, :math:\mathrm{mode} = 2 and :math:\mathrm{diag} contained a non-positive element.

(errno :math:18)
On entry, :math:\mathrm{maxfev} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{maxfev} > 0.

**Warns**
**NagAlgorithmicWarning**
(errno :math:2)
There have been at least :math:\mathrm{maxfev} calls to :math:\mathrm{fcnval} and :math:\mathrm{fcnjac} combined: :math:\mathrm{maxfev} = \langle\mathit{\boldsymbol{value}}\rangle. Consider restarting the calculation from the final point held in :math:\mathrm{x}.

(errno :math:3)
No further improvement in the solution is possible. :math:\mathrm{xtol} is too small: :math:\mathrm{xtol} = \langle\mathit{\boldsymbol{value}}\rangle.

(errno :math:4)
The iteration is not making good progress, as measured by the improvement from the last :math:\langle\mathit{\boldsymbol{value}}\rangle Jacobian evaluations. This failure exit may indicate that the system does not have a zero, or that the solution is very close to the origin (see Accuracy <https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05rcf.html#accuracy>__). Otherwise, rerunning sys_deriv_expert from a different starting point may avoid the region of difficulty.

(errno :math:5)
The iteration is not making good progress, as measured by the improvement from the last :math:\langle\mathit{\boldsymbol{value}}\rangle iterations. This failure exit may indicate that the system does not have a zero, or that the solution is very close to the origin (see Accuracy <https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05rcf.html#accuracy>__). Otherwise, rerunning sys_deriv_expert from a different starting point may avoid the region of difficulty.

**NagCallbackTerminateWarning**
(errno :math:6)
Termination requested in :math:\mathrm{fcnval}, :math:\mathrm{fcnjac} or :math:\mathrm{fcnmon}.

.. _c05rc-py2-py-notes:

**Notes**
The system of equations is defined as:

.. math::
f_i\left(x_1, x_2, \ldots, x_n\right) = 0\text{, }i = 1,2,\ldots,n\text{.}

sys_deriv_expert is based on the MINPACK routine HYBRJ (see Moré et al. (1980)).
It chooses the correction at each step as a convex combination of the Newton and scaled gradient directions.
The Jacobian is updated by the rank-1 method of Broyden.
At the starting point, the Jacobian is requested, but it is not asked for again until the rank-1 method fails to produce satisfactory progress.
For more details see Powell (1970).

.. _c05rc-py2-py-references:

**References**
Moré, J J, Garbow, B S and Hillstrom, K E, 1980, User guide for MINPACK-1, Technical Report ANL-80-74, Argonne National Laboratory

Powell, M J D, 1970, A hybrid method for nonlinear algebraic equations, Numerical Methods for Nonlinear Algebraic Equations, (ed P Rabinowitz), Gordon and Breach
"""
raise NotImplementedError

[docs]def sys_deriv_rcomm(irevcm, x, fvec, fjac, mode, diag, r, qtf, comm, xtol=sqrt(machine.precision()), factor=100.0):
r"""
sys_deriv_rcomm is a comprehensive reverse communication function that finds a solution of a system of nonlinear equations by a modification of the Powell hybrid method.
You must provide the Jacobian.

.. _c05rd-py2-py-doc:

For full information please refer to the NAG Library document for c05rd

https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05rdf.html

.. _c05rd-py2-py-parameters:

**Parameters**
**irevcm** : int
On initial entry: must have the value :math:0.

**x** : float, ndarray, shape :math:\left(n\right), modified in place
On initial entry: an initial guess at the solution vector.

On intermediate exit: contains the current point.

On final exit: the final estimate of the solution vector.

**fvec** : float, ndarray, shape :math:\left(n\right), modified in place
On initial entry: need not be set.

On intermediate entry: if :math:\mathrm{irevcm}\neq 2, :math:\mathrm{fvec} must not be changed.

If :math:\mathrm{irevcm} = 2, :math:\mathrm{fvec} must be set to the values of the functions computed at the current point :math:\mathrm{x}.

On final exit: the function values at the final point, :math:\mathrm{x}.

**fjac** : float, ndarray, shape :math:\left(n, n\right), modified in place
On initial entry: need not be set.

On intermediate entry: if :math:\mathrm{irevcm}\neq 3, :math:\mathrm{fjac} must not be changed.

If :math:\mathrm{irevcm} = 3, :math:\mathrm{fjac}[\textit{i}-1,\textit{j}-1] must contain the value of :math:\frac{{\partial f_{\textit{i}}}}{{\partial x_{\textit{j}}}} at the point :math:x, for :math:\textit{j} = 1,2,\ldots,n, for :math:\textit{i} = 1,2,\ldots,n.

On final exit: the orthogonal matrix :math:Q produced by the :math:QR factorization of the final approximate Jacobian.

**mode** : int
On initial entry: indicates whether or not you have provided scaling factors in :math:\mathrm{diag}.

If :math:\mathrm{mode} = 2, the scaling must have been supplied in :math:\mathrm{diag}.

Otherwise, if :math:\mathrm{mode} = 1, the variables will be scaled internally.

**diag** : float, ndarray, shape :math:\left(n\right), modified in place
On initial entry: if :math:\mathrm{mode} = 2, :math:\mathrm{diag} must contain multiplicative scale factors for the variables.

If :math:\mathrm{mode} = 1, :math:\mathrm{diag} need not be set.

On intermediate exit: :math:\mathrm{diag} must not be changed.

On final exit: the scale factors actually used (computed internally if :math:\mathrm{mode} = 1).

**r** : float, ndarray, shape :math:\left(n\times \left(n+1\right)/2\right), modified in place
On initial entry: need not be set.

On intermediate exit: must not be changed.

On final exit: the upper triangular matrix :math:R produced by the :math:QR factorization of the final approximate Jacobian, stored row-wise.

**qtf** : float, ndarray, shape :math:\left(n\right), modified in place
On initial entry: need not be set.

On intermediate exit: must not be changed.

On final exit: the vector :math:Q^\mathrm{T}f.

**comm** : dict, communication object, modified in place
Communication structure.

On initial entry: need not be set.

**xtol** : float, optional
On initial entry: the accuracy in :math:\mathrm{x} to which the solution is required.

**factor** : float, optional
On initial entry: a quantity to be used in determining the initial step bound. In most cases, :math:\mathrm{factor} should lie between :math:0.1 and :math:100.0. (The step bound is :math:\mathrm{factor}\times \left\lVert \mathrm{diag}\times \mathrm{x}\right\rVert_2 if this is nonzero; otherwise the bound is :math:\mathrm{factor}.)

**Returns**
**irevcm** : int
On intermediate exit: specifies what action you must take before re-entering sys_deriv_rcomm with :math:\mathrm{irevcm} set to this value. The value of :math:\mathrm{irevcm} should be interpreted as follows:

:math:\mathrm{irevcm} = 1

Indicates the start of a new iteration. No action is required by you, but :math:\mathrm{x} and :math:\mathrm{fvec} are available for printing.

:math:\mathrm{irevcm} = 2

Indicates that before re-entry to sys_deriv_rcomm, :math:\mathrm{fvec} must contain the function values :math:f_i\left(x\right).

:math:\mathrm{irevcm} = 3

Indicates that before re-entry to sys_deriv_rcomm, :math:\mathrm{fjac}[\textit{i}-1,\textit{j}-1] must contain the value of :math:\frac{{\partial f_{\textit{i}}}}{{\partial x_{\textit{j}}}} at the point :math:x, for :math:\textit{j} = 1,2,\ldots,n, for :math:\textit{i} = 1,2,\ldots,n.

On final exit: :math:\mathrm{irevcm} = 0 and the algorithm has terminated.

.. _c05rd-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:2)
On entry, :math:\mathrm{irevcm} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{irevcm} = 0, :math:1, :math:2 or :math:3.

(errno :math:11)
On entry, :math:n = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:n > 0.

(errno :math:12)
On entry, :math:\mathrm{xtol} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{xtol}\geq 0.0.

(errno :math:13)
On entry, :math:\mathrm{mode} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{mode} = 1 or :math:2.

(errno :math:14)
On entry, :math:\mathrm{factor} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{factor} > 0.0.

(errno :math:15)
On entry, :math:\mathrm{mode} = 2 and :math:\mathrm{diag} contained a non-positive element.

**Warns**
**NagAlgorithmicWarning**
(errno :math:3)
No further improvement in the solution is possible. :math:\mathrm{xtol} is too small: :math:\mathrm{xtol} = \langle\mathit{\boldsymbol{value}}\rangle.

(errno :math:4)
The iteration is not making good progress, as measured by the improvement from the last :math:\langle\mathit{\boldsymbol{value}}\rangle Jacobian evaluations. This failure exit may indicate that the system does not have a zero, or that the solution is very close to the origin (see Accuracy <https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05rdf.html#accuracy>__). Otherwise, rerunning sys_deriv_rcomm from a different starting point may avoid the region of difficulty.

(errno :math:5)
The iteration is not making good progress, as measured by the improvement from the last :math:\langle\mathit{\boldsymbol{value}}\rangle iterations. This failure exit may indicate that the system does not have a zero, or that the solution is very close to the origin (see Accuracy <https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05rdf.html#accuracy>__). Otherwise, rerunning sys_deriv_rcomm from a different starting point may avoid the region of difficulty.

.. _c05rd-py2-py-notes:

**Notes**
The system of equations is defined as:

.. math::
f_i\left(x_1, x_2, \ldots, x_n\right) = 0\text{, }\quad i = 1,2,\ldots,n\text{.}

sys_deriv_rcomm is based on the MINPACK routine HYBRJ (see Moré et al. (1980)).
It chooses the correction at each step as a convex combination of the Newton and scaled gradient directions.
The Jacobian is updated by the rank-1 method of Broyden.
For more details see Powell (1970).

.. _c05rd-py2-py-references:

**References**
Moré, J J, Garbow, B S and Hillstrom, K E, 1980, User guide for MINPACK-1, Technical Report ANL-80-74, Argonne National Laboratory

Powell, M J D, 1970, A hybrid method for nonlinear algebraic equations, Numerical Methods for Nonlinear Algebraic Equations, (ed P Rabinowitz), Gordon and Breach
"""
raise NotImplementedError

[docs]def sys_deriv_check(mode, x, fvec, fjac, fvecp):
r"""
sys_deriv_check checks the user-supplied gradients of a set of nonlinear functions in several variables, for consistency with the functions themselves.
The function must be called twice.

.. _c05zd-py2-py-doc:

For full information please refer to the NAG Library document for c05zd

https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05zdf.html

.. _c05zd-py2-py-parameters:

**Parameters**
**mode** : int
The value :math:1 on the first call and the value :math:2 on the second call of sys_deriv_check.

**x** : float, array-like, shape :math:\left(n\right)
The components of a point :math:x, at which the consistency check is to be made. (See Accuracy <https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05zdf.html#accuracy>__.)

**fvec** : float, array-like, shape :math:\left(m\right)
If :math:\mathrm{mode} = 2, :math:\mathrm{fvec} must contain the value of the functions evaluated at :math:x. If :math:\mathrm{mode} = 1, :math:\mathrm{fvec} is not referenced.

**fjac** : float, array-like, shape :math:\left(m, n\right)
If :math:\mathrm{mode} = 2, :math:\mathrm{fjac} must contain the value of :math:\frac{{\partial f_i}}{{\partial x_j}} at the point :math:x, for :math:\textit{j} = 1,2,\ldots,n, for :math:\textit{i} = 1,2,\ldots,m. If :math:\mathrm{mode} = 1, :math:\mathrm{fjac} is not referenced.

**fvecp** : float, array-like, shape :math:\left(m\right)
If :math:\mathrm{mode} = 2, :math:\mathrm{fvecp} must contain the value of the functions evaluated at :math:\mathrm{xp} (as output by a preceding call to sys_deriv_check with :math:\mathrm{mode} = 1). If :math:\mathrm{mode} = 1, :math:\mathrm{fvecp} is not referenced.

**Returns**
**xp** : float, ndarray, shape :math:\left(n\right)
If :math:\mathrm{mode} = 1, :math:\mathrm{xp} is set to a point neighbouring :math:\mathrm{x}. If :math:\mathrm{mode} = 2, :math:\mathrm{xp} is undefined.

**err** : float, ndarray, shape :math:\left(m\right)
If :math:\mathrm{mode} = 2, :math:\mathrm{err} contains measures of correctness of the respective gradients. If :math:\mathrm{mode} = 1, :math:\mathrm{err} is undefined. If there is no loss of significance (see Accuracy <https://www.nag.com/numeric/nl/nagdoc_29.2/flhtml/c05/c05zdf.html#accuracy>__), if :math:\mathrm{err}[i-1] is :math:1.0 the :math:i\ th user-supplied gradient :math:\frac{{\partial f_i}}{{\partial x_j}}, for :math:\textit{j} = 1,2,\ldots,n is correct, whilst if :math:\mathrm{err}[i-1] is :math:0.0 the :math:i\ th gradient is incorrect. For values of :math:\mathrm{err}[i-1] between :math:0.0 and :math:1.0 the categorisation is less certain. In general, a value of :math:\mathrm{err}[i-1] > 0.5 indicates that the :math:i\ th gradient is probably correct.

.. _c05zd-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:1)
On entry, :math:\mathrm{mode} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{mode} = 1 or :math:2.

(errno :math:2)
On entry, :math:m = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:m \geq 1.

(errno :math:3)
On entry, :math:n = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:n \geq 1.

.. _c05zd-py2-py-notes:

**Notes**
sys_deriv_check is based on the MINPACK routine CHKDER (see Moré et al. (1980)).
It checks the :math:i\ th gradient for consistency with the :math:i\ th function by computing a forward-difference approximation along a suitably chosen direction and comparing this approximation with the user-supplied gradient along the same direction.
The principal characteristic of sys_deriv_check is its invariance under changes in scale of the variables or functions.

.. _c05zd-py2-py-references:

**References**
Moré, J J, Garbow, B S and Hillstrom, K E, 1980, User guide for MINPACK-1, Technical Report ANL-80-74, Argonne National Laboratory
"""
raise NotImplementedError