# Source code for naginterfaces.library.numdiff

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

numdiff - Numerical Differentiation

This module is concerned with calculating approximations to derivatives of a function :math:f.

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

Generates abscissae for :meth:rcomm: :meth:sample

**Numerical derivatives**

direct communication: :meth:fwd

reverse communication: :meth:rcomm

For full information please refer to the NAG Library document

https://www.nag.com/numeric/nl/nagdoc_28.4/flhtml/d04/d04intro.html
"""

[docs]def fwd(xval, nder, hbase, f, data=None):
r"""
fwd calculates a set of derivatives (up to order :math:14) of a function of one real variable at a point, together with a corresponding set of error estimates, using an extension of the Neville algorithm.

.. _d04aa-py2-py-doc:

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

https://www.nag.com/numeric/nl/nagdoc_28.4/flhtml/d04/d04aaf.html

.. _d04aa-py2-py-parameters:

**Parameters**
**xval** : float
The point at which the derivatives are required, :math:x_0.

**nder** : int
Must be set so that its absolute value is the highest order derivative required.

:math:\mathrm{nder} > 0

All derivatives up to order :math:\mathrm{min}\left(\mathrm{nder}, 14\right) are calculated.

:math:\mathrm{nder} < 0 and :math:\mathrm{nder} is even

Only even order derivatives up to order :math:\mathrm{min}\left({-\mathrm{nder}}, 14\right) are calculated.

:math:\mathrm{nder} < 0 and :math:\mathrm{nder} is odd

Only odd order derivatives up to order :math:\mathrm{min}\left({-\mathrm{nder}}, 13\right) are calculated.

**hbase** : float
The initial step length which may be positive or negative. For advice on the choice of :math:\mathrm{hbase} see Further Comments <https://www.nag.com/numeric/nl/nagdoc_28.4/flhtml/d04/d04aaf.html#fcomments>__.

**f** : callable retval = f(x, data=None)
:math:\mathrm{f} must evaluate the function :math:f\left(x\right) at a specified point.

**Parameters**
**x** : float
The value of the argument :math:x.

If you have equally spaced tabular data, the following information may be useful:

(i) in any call of fwd the only values of :math:x for which :math:f\left(x\right) will be required are :math:x = \mathrm{xval} and :math:x = \mathrm{xval}\pm \left(2\textit{j}-1\right)\mathrm{hbase}, for :math:\textit{j} = 1,2,\ldots,10; and

(#) :math:f\left(x_0\right) is always computed, but it is disregarded when only odd order derivatives are required.

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

**Returns**
**retval** : float
The value of :math:f\left(x\right) at the specified point.

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

**Returns**
**der** : float, ndarray, shape :math:\left(14\right)
:math:\mathrm{der}[j-1] contains an approximation to the :math:j\ th derivative of :math:f\left(x\right) at :math:x = \mathrm{xval}, so long as the :math:j\ th derivative is one of those requested by you when specifying :math:\mathrm{nder}. For other values of :math:j, :math:\mathrm{der}[j-1] is unused.

**erest** : float, ndarray, shape :math:\left(14\right)
An estimate of the absolute error in the corresponding result :math:\mathrm{der}[j-1] so long as the :math:j\ th derivative is one of those requested by you when specifying :math:\mathrm{nder}. The sign of :math:\mathrm{erest}[j-1] is positive unless the result :math:\mathrm{der}[j-1] is questionable. It is set negative when :math:\left\lvert \mathrm{der}[j-1]\right\rvert < \left\lvert \mathrm{erest}[j-1]\right\rvert or when for some other reason there is doubt about the validity of the result :math:\mathrm{der}[j-1] (see :ref:Exceptions <d04aa-py2-py-errors>). For other values of :math:j, :math:\mathrm{erest}[j-1] is unused.

.. _d04aa-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:1)
On entry, :math:\mathrm{nder} = 0.

Constraint: :math:\mathrm{nder} \neq 0.

(errno :math:1)
On entry, :math:\mathrm{hbase} = 0.0.

Constraint: :math:\mathrm{hbase}\neq 0.0.

.. _d04aa-py2-py-notes:

**Notes**
fwd provides a set of approximations:

.. math::

to the derivatives:

.. math::

of a real valued function :math:f\left(x\right) at a real abscissa :math:x_0, together with a set of error estimates:

.. math::

which hopefully satisfy:

.. math::
\left\lvert \mathrm{der}[j-1]-f^{\left(j\right)}\left(x_0\right)\right\rvert < \mathrm{erest}[j-1]\text{, }\quad j = 1,2,\ldots,n\text{.}

You must provide the value of :math:x_0, a value of :math:n (which is reduced to :math:14 should it exceed :math:14), a function which evaluates :math:f\left(x\right) for all real :math:x, and a step length :math:h.
The results :math:\mathrm{der}[j-1] and :math:\mathrm{erest}[j-1] are based on :math:21 function values:

.. math::
f\left(x_0\right),f\left(x_0\pm \left(2i-1\right)h\right)\text{, }\quad i = 1,2,\ldots,10\text{.}

Internally fwd calculates the odd order derivatives and the even order derivatives separately.
There is an option you can use for restricting the calculation to only odd (or even) order derivatives.
For each derivative the function employs an extension of the Neville Algorithm (see Lyness and Moler (1969)) to obtain a selection of approximations.

For example, for odd derivatives, based on :math:20 function values, fwd calculates a set of numbers:

.. math::

each of which is an approximation to :math:f^{\left(2s+1\right)}\left(x_0\right)/\left(2s+1\right)!.
A specific approximation :math:T_{{\textit{k},p,s}} is of polynomial degree :math:2p+2 and is based on polynomial interpolation using function values :math:f\left(x_0\pm \left(2i-1\right)h\right), for :math:\textit{k} = \textit{k},\ldots,\textit{k}+p.
In the absence of round-off error, the better approximations would be associated with the larger values of :math:p and of :math:k.
However, round-off error in function values has an increasingly contaminating effect for successively larger values of :math:p.
This function proceeds to make a judicious choice between all the approximations in the following way.

For a specified value of :math:s, let:

.. math::
R_p = U_p-L_p\text{, }\quad p = s,s+1,\ldots,6

where :math:U_p = \mathrm{max}_{\textit{k}}\left(T_{{\textit{k},p,s}}\right) and :math:L_p = \mathrm{min}_{\textit{k}}\left(T_{{\textit{k},p,s}}\right), for :math:\textit{k} = 0,1,\ldots,9-p, and let :math:\bar{\textit{p}} be such that :math:R_{{\bar{\textit{p}}}} = \mathrm{min}_{\textit{p}}\left(R_{\textit{p}}\right), for :math:\textit{p} = s,\ldots,6.

The function returns:

.. math::
\mathrm{der}[2s] = \frac{1}{{8-\bar{p}}}\times \left\{\sum_{{k = 0}}^{{9-\bar{p}}}T_{{k,\bar{p},s}}-U_{{\bar{p}}}-L_{{\bar{p}}}\right\}\left(2s+1\right)!

and

.. math::
\mathrm{erest}[2s] = R_{\bar{p}}\times \left(2s+1\right)!\times K_{{2s+1}}

where :math:K_j is a safety factor which has been assigned the values:

.. rst-class:: nag-rules-none

+------------------+-----------------+
|:math:K_j = 1,  |:math:j\leq 9  |
+------------------+-----------------+
|:math:K_j = 1.5,|:math:j = 10,11|
+------------------+-----------------+
|:math:K_j = 2   |:math:j\geq 12,|
+------------------+-----------------+

on the basis of performance statistics.

The even order derivatives are calculated in a precisely analogous manner.

.. _d04aa-py2-py-references:

**References**
Lyness, J N and Moler, C B, 1966, van der Monde systems and numerical differentiation, Numer. Math. (8), 458--464

Lyness, J N and Moler, C B, 1969, Generalised Romberg methods for integrals of derivatives, Numer. Math. (14), 1--14
"""
raise NotImplementedError

[docs]def rcomm(xval, fval):
r"""
rcomm calculates a set of derivatives (up to order :math:14) of a function at a point with respect to a single variable.
A corresponding set of error estimates is also returned.
Derivatives are calculated using an extension of the Neville algorithm.
This function differs from :meth:fwd, in that the abscissae and corresponding function values must be calculated before this function is called.
The abscissae may be generated using :meth:sample.

.. _d04ba-py2-py-doc:

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

https://www.nag.com/numeric/nl/nagdoc_28.4/flhtml/d04/d04baf.html

.. _d04ba-py2-py-parameters:

**Parameters**
**xval** : float, array-like, shape :math:\left(21\right)
The abscissae at which the function has been evaluated, as described in :ref:Notes <d04ba-py2-py-notes>. These can be generated by calling :meth:sample. The order of the abscissae is irrelevant.

**fval** : float, array-like, shape :math:\left(21\right)
:math:\mathrm{fval}[\textit{j}-1] must contain the function value at :math:\mathrm{xval}[\textit{j}-1], for :math:\textit{j} = 1,2,\ldots,21.

**Returns**
**der** : float, ndarray, shape :math:\left(14\right)
The :math:14 derivative estimates.

**erest** : float, ndarray, shape :math:\left(14\right)
The :math:14 error estimates for the derivatives.

.. _d04ba-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:1)
The derived :math:h is below tolerance.

Derived :math:h > \langle\mathit{\boldsymbol{value}}\rangle is required. Derived :math:h = \langle\mathit{\boldsymbol{value}}\rangle.

(errno :math:1)
On entry, the values of :math:\mathrm{xval} are not correctly spaced.

Derived :math:h = \langle\mathit{\boldsymbol{value}}\rangle.

.. _d04ba-py2-py-notes:

**Notes**
rcomm provides a set of approximations:

.. math::

to the derivatives:

.. math::

of a real valued function :math:f\left(x\right) at a real abscissa :math:x_0, together with a set of error estimates:

.. math::

which hopefully satisfy:

.. math::
\left\lvert \mathrm{der}[j-1]-f^{\left(j\right)}\left(x_0\right)\right\rvert < \mathrm{erest}[j-1]\text{, }\quad j = 1,2,\ldots,14\text{.}

The results :math:\mathrm{der}[j-1] and :math:\mathrm{erest}[j-1] are based on :math:21 function values:

.. math::
f\left(x_0\right),f\left(x_0\pm \left(2i-1\right)h\right)\text{, }\quad i = 1,2,\ldots,10\text{.}

The abscissae :math:x and the corresponding function values :math:f\left(x\right) should be passed into rcomm as the vectors :math:\mathrm{xval} and :math:\mathrm{fval} respectively.
The step size :math:h is derived from the abscissae in :math:\mathrm{xval}.
See Further Comments <https://www.nag.com/numeric/nl/nagdoc_28.4/flhtml/d04/d04baf.html#fcomments>__ for a discussion of how the derived value of :math:h may affect the results of rcomm.
The order in which the abscissae and function values are stored in :math:\mathrm{xval} and :math:\mathrm{fval} is irrelevant, provided that the function value at any given index corresponds to the value of the abscissa at the same index.
Abscissae may be automatically generated using :meth:sample if desired.
For each derivative rcomm employs an extension of the Neville Algorithm (see Lyness and Moler (1969)) to obtain a selection of approximations.

For example, for odd derivatives, this function calculates a set of numbers:

.. math::

each of which is an approximation to :math:f^{\left(2s+1\right)}\left(x_0\right)/\left(2s+1\right)!.
A specific approximation :math:T_{{\textit{k},p,s}} is of polynomial degree :math:2p+2 and is based on polynomial interpolation using function values :math:f\left(x_0\pm \left(2i-1\right)h\right), for :math:\textit{k} = \textit{k},\ldots,\textit{k}+p.
In the absence of round-off error, the better approximations would be associated with the larger values of :math:p and of :math:k.
However, round-off error in function values has an increasingly contaminating effect for successively larger values of :math:p.
This function proceeds to make a judicious choice between all the approximations in the following way.

For a specified value of :math:s, let:

.. math::
R_p = U_p-L_p\text{, }\quad p = s,s+1,\ldots,6

where :math:U_p = \mathrm{max}_{\textit{k}}\left(T_{{\textit{k},p,s}}\right) and :math:L_p = \mathrm{min}_{\textit{k}}\left(T_{{\textit{k},p,s}}\right), for :math:\textit{k} = 0,1,\ldots,9-p, and let :math:\bar{p} be such that :math:R_{{\bar{\textit{p}}}} = \mathrm{min}_{\textit{p}}\left(R_{\textit{p}}\right), for :math:\textit{p} = s,\ldots,6.

This function returns:

.. math::
\mathrm{der}[2s] = \frac{1}{{8-\bar{p}}}\times \left\{\sum_{{k = 0}}^{{9-\bar{p}}}T_{{k,\bar{p},s}}-U_{{\bar{p}}}-L_{{\bar{p}}}\right\}\left(2s+1\right)!

and

.. math::
\mathrm{erest}[2s] = R_{\bar{p}}\times \left(2s+1\right)!\times K_{{2s+1}}

where :math:K_j is a safety factor which has been assigned the values:

.. rst-class:: nag-rules-none

+------------------+-----------------+
|:math:K_j = 1,  |:math:j\leq 9  |
+------------------+-----------------+
|:math:K_j = 1.5,|:math:j = 10,11|
+------------------+-----------------+
|:math:K_j = 2   |:math:j\geq 12,|
+------------------+-----------------+

on the basis of performance statistics.

The even order derivatives are calculated in a precisely analogous manner.

.. _d04ba-py2-py-references:

**References**
Lyness, J N and Moler, C B, 1969, Generalised Romberg methods for integrals of derivatives, Numer. Math. (14), 1--14
"""
raise NotImplementedError

[docs]def sample(x_0, hbase):
r"""
sample generates abscissae about a target abscissa :math:x_0 for use in a subsequent call to :meth:rcomm.

.. _d04bb-py2-py-doc:

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

https://www.nag.com/numeric/nl/nagdoc_28.4/flhtml/d04/d04bbf.html

.. _d04bb-py2-py-parameters:

**Parameters**
**x_0** : float
The abscissa :math:x_0 at which derivatives are required.

**hbase** : float
The chosen step size :math:h. If :math:h < 10\epsilon, where :math:\epsilon = \texttt{machine.precision}\left(\right), the default :math:h = \epsilon^{\left(1/4\right)} will be used.

**Returns**
**xval** : float, ndarray, shape :math:\left(21\right)
The abscissae for passing to :meth:rcomm.

.. _d04bb-py2-py-notes:

**Notes**
sample may be used to generate the necessary abscissae about a target abscissa :math:x_0 for the calculation of derivatives using :meth:rcomm.

For a given :math:x_0 and :math:h, the abscissae correspond to the set :math:\left\{x_0, {x_0\pm \left(2\textit{j}-1\right)h}\right\}, for :math:\textit{j} = 1,2,\ldots,10.
These :math:21 points will be returned in ascending order in :math:\mathrm{xval}.
In particular, :math:\mathrm{xval}[10] will be equal to :math:x_0.

.. _d04bb-py2-py-references:

**References**
Lyness, J N and Moler, C B, 1969, Generalised Romberg methods for integrals of derivatives, Numer. Math. (14), 1--14
"""
raise NotImplementedError