# -*- coding: utf-8 -*-
r"""
Module Summary
--------------
Interfaces for the NAG Mark 30.1 `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::
f\left(w\right) = we^w\quad \text{ for }\quad w \in \mathbb{C}\text{;}
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://support.nag.com/numeric/nl/nagdoc_30.1/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[0]
... 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[0] = (2.*x[1] - 1.)/(3. - 2.*x[0])
... 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,
)
See Also
--------
``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://support.nag.com/numeric/nl/nagdoc_30.1/flhtml/c05/c05intro.html
"""
# NAG Copyright 2017-2024.
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://support.nag.com/numeric/nl/nagdoc_30.1/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://support.nag.com/numeric/nl/nagdoc_30.1/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
See Also
--------
: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://support.nag.com/numeric/nl/nagdoc_30.1/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://support.nag.com/numeric/nl/nagdoc_30.1/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://support.nag.com/numeric/nl/nagdoc_30.1/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://support.nag.com/numeric/nl/nagdoc_30.1/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::
f\left(x\right)-\theta_rf\left(x_0\right)\text{, }\quad r = 0,1,\ldots,m
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://support.nag.com/numeric/nl/nagdoc_30.1/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://support.nag.com/numeric/nl/nagdoc_30.1/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://support.nag.com/numeric/nl/nagdoc_30.1/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}[0]` 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}[0]` need not be set.
:math:`\mathrm{ind} = -1`
:math:`\mathrm{fx}` and :math:`\mathrm{c}[0]` 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://support.nag.com/numeric/nl/nagdoc_30.1/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::
f\left(w\right) = we^w\quad \text{ for }\quad w \in C\text{.}
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://support.nag.com/numeric/nl/nagdoc_30.1/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::
f\left(w\right) = we^w\quad \text{ for }\quad w \in C\text{.}
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://support.nag.com/numeric/nl/nagdoc_30.1/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://support.nag.com/numeric/nl/nagdoc_30.1/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://support.nag.com/numeric/nl/nagdoc_30.1/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://support.nag.com/numeric/nl/nagdoc_30.1/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://support.nag.com/numeric/nl/nagdoc_30.1/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`
No calls are made.
: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://support.nag.com/numeric/nl/nagdoc_30.1/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
See Also
--------
: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://support.nag.com/numeric/nl/nagdoc_30.1/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://support.nag.com/numeric/nl/nagdoc_30.1/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://support.nag.com/numeric/nl/nagdoc_30.1/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://support.nag.com/numeric/nl/nagdoc_30.1/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://support.nag.com/numeric/nl/nagdoc_30.1/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`
No calls are made.
: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://support.nag.com/numeric/nl/nagdoc_30.1/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://support.nag.com/numeric/nl/nagdoc_30.1/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://support.nag.com/numeric/nl/nagdoc_30.1/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://support.nag.com/numeric/nl/nagdoc_30.1/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://support.nag.com/numeric/nl/nagdoc_30.1/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://support.nag.com/numeric/nl/nagdoc_30.1/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://support.nag.com/numeric/nl/nagdoc_30.1/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://support.nag.com/numeric/nl/nagdoc_30.1/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
