# -*- coding: utf-8 -*-
r"""
Module Summary
--------------
Interfaces for the NAG Mark 30.3 `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://support.nag.com/numeric/nl/nagdoc_30.3/flhtml/d04/d04intro.html
"""
# NAG Copyright 2017-2024.
[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://support.nag.com/numeric/nl/nagdoc_30.3/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://support.nag.com/numeric/nl/nagdoc_30.3/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::
\mathrm{der}[j-1]\text{, }\quad j = 1,2,\ldots,n
to the derivatives:
.. math::
f^{\left(j\right)}\left(x_0\right)\text{, }\quad j = 1,2,\ldots,n
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::
\mathrm{erest}[j-1]\text{, }\quad j = 1,2,\ldots,n
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::
T_{{k,p,s}}\text{, }\quad p = s,s+1,\ldots,6\text{, }\quad k = 0,1,\ldots,9-p
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://support.nag.com/numeric/nl/nagdoc_30.3/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::
\mathrm{der}[j-1]\text{, }\quad j = 1,2,\ldots,14
to the derivatives:
.. math::
f^{\left(j\right)}\left(x_0\right)\text{, }\quad j = 1,2,\ldots,14
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::
\mathrm{erest}[j-1]\text{, }\quad j = 1,2,\ldots,14
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://support.nag.com/numeric/nl/nagdoc_30.3/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::
T_{{k,p,s}}\text{, }\quad p = s,s+1,\ldots,6\text{, }\quad k = 0,1,\ldots,9-p
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://support.nag.com/numeric/nl/nagdoc_30.3/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
