# -*- coding: utf-8 -*-
r"""
Module Summary
--------------
Interfaces for the NAG Mark 30.2 `inteq` Chapter.
``inteq`` - Integral Equations
This module is concerned with the numerical solution of integral equations.
Provision will be made for most of the standard types of equation (see `Background to the Problems <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05intro.html#background>`__).
The following are, however, specifically excluded:
(a) Equations arising in the solution of partial differential equations by integral equation methods. In cases where the prime purpose of an algorithm is the solution of a partial differential equation it will normally be included in submodule :mod:`~naginterfaces.library.pde`.
(#) Calculation of inverse integral transforms. This problem falls within the scope of submodule :mod:`~naginterfaces.library.sum`.
Functionality Index
-------------------
**Fredholm equation of second kind**
linear
nonsingular discontinuous or 'split' kernel: :meth:`fredholm2_split`
nonsingular smooth kernel: :meth:`fredholm2_smooth`
**Volterra equation of first kind**
nonlinear
weakly-singular
convolution equation (Abel): :meth:`abel1_weak`
**Volterra equation of second kind**
nonlinear
nonsingular
convolution equation: :meth:`volterra2`
weakly-singular
convolution equation (Abel): :meth:`abel2_weak`
**Weight generating functions**
weights for general solution of Volterra equations: :meth:`volterra_weights`
weights for general solution of Volterra equations with weakly-singular kernel: :meth:`abel_weak_weights`
For full information please refer to the NAG Library document
https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05intro.html
"""
# NAG Copyright 2017-2024.
from __future__ import absolute_import
from math import sqrt
from ..library import machine
[docs]def fredholm2_split(lamda, a, b, k1, k2, g, n, ind, data=None):
r"""
``fredholm2_split`` solves a linear, nonsingular Fredholm equation of the second kind with a split kernel.
.. _d05aa-py2-py-doc:
For full information please refer to the NAG Library document for d05aa
https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05aaf.html
.. _d05aa-py2-py-parameters:
**Parameters**
**lamda** : float
The value of the parameter :math:`\lambda` of the integral equation.
**a** : float
:math:`a`, the lower limit of integration.
**b** : float
:math:`b`, the upper limit of integration.
**k1** : callable retval = k1(x, s, data=None)
:math:`\mathrm{k1}` must evaluate the kernel :math:`k\left(x, s\right) = k_1\left(x, s\right)` of the integral equation for :math:`a\leq s\leq x`.
**Parameters**
**x** : float
The values of :math:`x` and :math:`s` at which :math:`k_1\left(x, s\right)` is to be evaluated.
**s** : float
The values of :math:`x` and :math:`s` at which :math:`k_1\left(x, s\right)` is to be evaluated.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**Returns**
**retval** : float
The value of the kernel :math:`k\left(x, s\right) = k_1\left(x, s\right)` evaluated at :math:`\mathrm{x}` and :math:`\mathrm{s}`.
**k2** : callable retval = k2(x, s, data=None)
:math:`\mathrm{k2}` must evaluate the kernel :math:`k\left(x, s\right) = k_2\left(x, s\right)` of the integral equation for :math:`x < s\leq b`.
**Parameters**
**x** : float
The values of :math:`x` and :math:`s` at which :math:`k_2\left(x, s\right)` is to be evaluated.
**s** : float
The values of :math:`x` and :math:`s` at which :math:`k_2\left(x, s\right)` is to be evaluated.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**Returns**
**retval** : float
The value of the kernel :math:`k\left(x, s\right) = k_2\left(x, s\right)` evaluated at :math:`\mathrm{x}` and :math:`\mathrm{s}`.
**g** : callable retval = g(x, data=None)
:math:`\mathrm{g}` must evaluate the function :math:`g\left(x\right)` for :math:`a\leq x\leq b`.
**Parameters**
**x** : float
The values of :math:`x` at which :math:`g\left(x\right)` is to be evaluated.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**Returns**
**retval** : float
The value of :math:`g\left(x\right)` evaluated at :math:`\mathrm{x}`.
**n** : int
The number of terms in the Chebyshev series required to approximate :math:`f\left(x\right)`.
**ind** : int
Determines the forms of the kernel, :math:`k\left(x, s\right)`, and the function :math:`g\left(x\right)`.
:math:`\mathrm{ind} = 0`
:math:`k\left(x, s\right)` is not centro-symmetric (or no account is to be taken of centro-symmetry).
:math:`\mathrm{ind} = 1`
:math:`k\left(x, s\right)` is centro-symmetric and :math:`g\left(x\right)` is odd.
:math:`\mathrm{ind} = 2`
:math:`k\left(x, s\right)` is centro-symmetric and :math:`g\left(x\right)` is even.
:math:`\mathrm{ind} = 3`
:math:`k\left(x, s\right)` is centro-symmetric but :math:`g\left(x\right)` is neither odd nor even.
**data** : arbitrary, optional
User-communication data for callback functions.
**Returns**
**f** : float, ndarray, shape :math:`\left(\mathrm{n}\right)`
The approximate values :math:`f_{\textit{i}}`, for :math:`\textit{i} = 1,2,\ldots,\mathrm{n}`, of :math:`f\left(x\right)` evaluated at the first :math:`\mathrm{n}` of :math:`m+1` Chebyshev points :math:`x_i`, (see :ref:`Notes <d05aa-py2-py-notes>`).
If :math:`\mathrm{ind} = 0` or :math:`3`, :math:`m = \mathrm{n}-1`.
If :math:`\mathrm{ind} = 1`, :math:`m = 2\times \mathrm{n}`.
If :math:`\mathrm{ind} = 2`, :math:`m = 2\times \mathrm{n}-1`.
**c** : float, ndarray, shape :math:`\left(\mathrm{n}\right)`
The coefficients :math:`c_{\textit{i}}`, for :math:`\textit{i} = 1,2,\ldots,\mathrm{n}`, of the Chebyshev series approximation to :math:`f\left(x\right)`.
If :math:`\mathrm{ind} = 1` this series contains polynomials of odd order only and if :math:`\mathrm{ind} = 2` the series contains even order polynomials only.
.. _d05aa-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{n}\geq 1`.
(`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{b} > \mathrm{a}`.
(`errno` :math:`2`)
A failure has occurred due to proximity of an eigenvalue.
.. _d05aa-py2-py-notes:
**Notes**
``fredholm2_split`` solves an integral equation of the form
.. math::
f\left(x\right)-\lambda \int_a^bk\left(x, s\right)f\left(s\right){ds} = g\left(x\right)
for :math:`a\leq x\leq b`, when the kernel :math:`k` is defined in two parts: :math:`k = k_1` for :math:`a\leq s\leq x` and :math:`k = k_2` for :math:`x < s\leq b`.
The method used is that of El--Gendi (1969) for which, it is important to note, each of the functions :math:`k_1` and :math:`k_2` must be defined, smooth and nonsingular, for all :math:`x` and :math:`s` in the interval :math:`\left[a, b\right]`.
An approximation to the solution :math:`f\left(x\right)` is found in the form of an :math:`n` term Chebyshev series :math:`{\sum^\prime}_{{i = 1}}^nc_iT_i\left(x\right)`, where :math:`{}^{\prime }` indicates that the first term is halved in the sum.
The coefficients :math:`c_{\textit{i}}`, for :math:`\textit{i} = 1,2,\ldots,n`, of this series are determined directly from approximate values :math:`f_{\textit{i}}`, for :math:`\textit{i} = 1,2,\ldots,n`, of the function :math:`f\left(x\right)` at the first :math:`n` of a set of :math:`m+1` Chebyshev points:
.. math::
x_i = \frac{1}{2}\left(a+b+\left(b-a\right)\cos\left[\left(i-1\right)\pi /m\right]\right)\text{, }\quad i = 1,2,\ldots,m+1\text{.}
The values :math:`f_i` are obtained by solving simultaneous linear algebraic equations formed by applying a quadrature formula (equivalent to the scheme of Clenshaw and Curtis (1960)) to the integral equation at the above points.
In general :math:`m = n-1`.
However, if the kernel :math:`k` is centro-symmetric in the interval :math:`\left[a, b\right]`, i.e., if :math:`k\left(x, s\right) = k\left({a+b-x}, {a+b-s}\right)`, then the function is designed to take advantage of this fact in the formation and solution of the algebraic equations.
In this case, symmetry in the function :math:`g\left(x\right)` implies symmetry in the function :math:`f\left(x\right)`.
In particular, if :math:`g\left(x\right)` is even about the mid-point of the range of integration, then so also is :math:`f\left(x\right)`, which may be approximated by an even Chebyshev series with :math:`m = 2n-1`.
Similarly, if :math:`g\left(x\right)` is odd about the mid-point then :math:`f\left(x\right)` may be approximated by an odd series with :math:`m = 2n`.
.. _d05aa-py2-py-references:
**References**
Clenshaw, C W and Curtis, A R, 1960, `A method for numerical integration on an automatic computer`, Numer. Math. (2), 197--205
El--Gendi, S E, 1969, `Chebyshev solution of differential, integral and integro-differential equations`, Comput. J. (12), 282--287
"""
raise NotImplementedError
[docs]def fredholm2_smooth(k, g, lamda, a, b, odorev, ev, n, data=None):
r"""
``fredholm2_smooth`` solves any linear nonsingular Fredholm integral equation of the second kind with a smooth kernel.
.. _d05ab-py2-py-doc:
For full information please refer to the NAG Library document for d05ab
https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05abf.html
.. _d05ab-py2-py-parameters:
**Parameters**
**k** : callable retval = k(x, s, data=None)
:math:`\mathrm{k}` must compute the value of the kernel :math:`k\left(x, s\right)` of the integral equation over the square :math:`a\leq x\leq b`, :math:`a\leq s\leq b`.
**Parameters**
**x** : float
The values of :math:`x` and :math:`s` at which :math:`k\left(x, s\right)` is to be calculated.
**s** : float
The values of :math:`x` and :math:`s` at which :math:`k\left(x, s\right)` is to be calculated.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**Returns**
**retval** : float
The value of the kernel :math:`k\left(x, s\right)` evaluated at :math:`\mathrm{x}` and :math:`\mathrm{s}`.
**g** : callable retval = g(x, data=None)
:math:`\mathrm{g}` must compute the value of the function :math:`g\left(x\right)` of the integral equation in the interval :math:`a\leq x\leq b`.
**Parameters**
**x** : float
The value of :math:`x` at which :math:`g\left(x\right)` is to be calculated.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**Returns**
**retval** : float
The value of :math:`g\left(x\right)` evaluated at :math:`\mathrm{x}`.
**lamda** : float
The value of the parameter :math:`\lambda` of the integral equation.
**a** : float
:math:`a`, the lower limit of integration.
**b** : float
:math:`b`, the upper limit of integration.
**odorev** : bool
Indicates whether it is known that the solution :math:`f\left(x\right)` is odd or even about the mid-point of the range of integration. If :math:`\mathrm{odorev}` is :math:`\mathbf{True}` then an odd or even solution is sought depending upon the value of :math:`\mathrm{ev}`.
**ev** : bool
Is ignored if :math:`\mathrm{odorev}` is :math:`\mathbf{False}`. Otherwise, if :math:`\mathrm{ev}` is :math:`\mathbf{True}`, an even solution is sought, whilst if :math:`\mathrm{ev}` is :math:`\mathbf{False}`, an odd solution is sought.
**n** : int
The number of terms in the Chebyshev series which approximates the solution :math:`f\left(x\right)`.
**data** : arbitrary, optional
User-communication data for callback functions.
**Returns**
**f** : float, ndarray, shape :math:`\left(\mathrm{n}\right)`
The approximate values :math:`f_{\textit{i}}`, for :math:`\textit{i} = 1,2,\ldots,\mathrm{n}`, of the function :math:`f\left(x\right)` at the first :math:`\mathrm{n}` of :math:`m+1` Chebyshev points (see :ref:`Notes <d05ab-py2-py-notes>`), where
.. rst-class:: nag-rules-none nag-align-left
+-------------------------+------------------------------------------------------------------------------------+
|:math:`m = 2\mathrm{n}-1`|if :math:`\mathrm{odorev} = \mathbf{True}` and :math:`\mathrm{ev} = \mathbf{True}`. |
+-------------------------+------------------------------------------------------------------------------------+
|:math:`m = 2\mathrm{n}` |if :math:`\mathrm{odorev} = \mathbf{True}` and :math:`\mathrm{ev} = \mathbf{False}`.|
+-------------------------+------------------------------------------------------------------------------------+
|:math:`m = \mathrm{n}-1` |if :math:`\mathrm{odorev} = \mathbf{False}`. |
+-------------------------+------------------------------------------------------------------------------------+
**c** : float, ndarray, shape :math:`\left(\mathrm{n}\right)`
The coefficients :math:`c_{\textit{i}}`, for :math:`\textit{i} = 1,2,\ldots,\mathrm{n}`, of the Chebyshev series approximation to :math:`f\left(x\right)`. When :math:`\mathrm{odorev}` is :math:`\mathbf{True}`, this series contains polynomials of even order only or of odd order only, according to :math:`\mathrm{ev}` being :math:`\mathbf{True}` or :math:`\mathbf{False}` respectively.
.. _d05ab-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\mathrm{n} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{n}\geq 1`.
(`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{b} > \mathrm{a}`.
(`errno` :math:`2`)
A failure has occurred due to proximity of an eigenvalue.
.. _d05ab-py2-py-notes:
**Notes**
``fredholm2_smooth`` uses the method of El--Gendi (1969) to solve an integral equation of the form
.. math::
f\left(x\right)-\lambda \int_a^bk\left(x, s\right)f\left(s\right){ds} = g\left(x\right)
for the function :math:`f\left(x\right)` in the range :math:`a\leq x\leq b`.
An approximation to the solution :math:`f\left(x\right)` is found in the form of an :math:`n` term Chebyshev series :math:`{\sum^\prime}_{{i = 1}}^nc_iT_i\left(x\right)`, where :math:`{}^{\prime }` indicates that the first term is halved in the sum.
The coefficients :math:`c_{\textit{i}}`, for :math:`\textit{i} = 1,2,\ldots,n`, of this series are determined directly from approximate values :math:`f_{\textit{i}}`, for :math:`\textit{i} = 1,2,\ldots,n`, of the function :math:`f\left(x\right)` at the first :math:`n` of a set of :math:`m+1` Chebyshev points
.. math::
x_i = \frac{1}{2}\left(a+b+\left(b-a\right)\times \cos\left[\left(i-1\right)\times \pi /m\right]\right)\text{, }\quad i = 1,2,\ldots,m+1\text{.}
The values :math:`f_i` are obtained by solving a set of simultaneous linear algebraic equations formed by applying a quadrature formula (equivalent to the scheme of Clenshaw and Curtis (1960)) to the integral equation at each of the above points.
In general :math:`m = n-1`.
However, advantage may be taken of any prior knowledge of the symmetry of :math:`f\left(x\right)`.
Thus if :math:`f\left(x\right)` is symmetric (i.e., even) about the mid-point of the range :math:`\left(a, b\right)`, it may be approximated by an even Chebyshev series with :math:`m = 2n-1`.
Similarly, if :math:`f\left(x\right)` is anti-symmetric (i.e., odd) about the mid-point of the range of integration, it may be approximated by an odd Chebyshev series with :math:`m = 2n`.
.. _d05ab-py2-py-references:
**References**
Clenshaw, C W and Curtis, A R, 1960, `A method for numerical integration on an automatic computer`, Numer. Math. (2), 197--205
El--Gendi, S E, 1969, `Chebyshev solution of differential, integral and integro-differential equations`, Comput. J. (12), 282--287
"""
raise NotImplementedError
[docs]def volterra2(ck, cg, cf, method, iorder, alim, tlim, nmesh, tol, thresh, lwk, data=None):
r"""
``volterra2`` computes the solution of a nonlinear convolution Volterra integral equation of the second kind using a reducible linear multi-step method.
.. _d05ba-py2-py-doc:
For full information please refer to the NAG Library document for d05ba
https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05baf.html
.. _d05ba-py2-py-parameters:
**Parameters**
**ck** : callable retval = ck(t, data=None)
:math:`\mathrm{ck}` must evaluate the kernel :math:`k\left(t\right)` of the integral equation `(1) <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05baf.html#eqn1>`__.
**Parameters**
**t** : float
:math:`t`, the value of the independent variable.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**Returns**
**retval** : float
The value of :math:`k` at the specified point.
**cg** : callable retval = cg(s, y, data=None)
:math:`\mathrm{cg}` must evaluate the function :math:`g\left(s, {y\left(s\right)}\right)` in `(1) <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05baf.html#eqn1>`__.
**Parameters**
**s** : float
:math:`s`, the value of the independent variable.
**y** : float
The value of the solution :math:`y` at the point :math:`\mathrm{s}`.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**Returns**
**retval** : float
The value of :math:`g` at the specified points.
**cf** : callable retval = cf(t, data=None)
:math:`\mathrm{cf}` must evaluate the function :math:`f\left(t\right)` in `(1) <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05baf.html#eqn1>`__.
**Parameters**
**t** : float
:math:`t`, the value of the independent variable.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**Returns**
**retval** : float
The value of :math:`f` at the specified point.
**method** : str, length 1
The type of method which you wish to employ.
:math:`\mathrm{method} = \texttt{'A'}`
For Adams' type formulae.
:math:`\mathrm{method} = \texttt{'B'}`
For backward differentiation formulae.
**iorder** : int
The order of the method to be used.
**alim** : float
:math:`a`, the lower limit of the integration interval.
**tlim** : float
The final point of the integration interval, :math:`T`.
**nmesh** : int
The number of equidistant points at which the solution is sought.
**tol** : float
The relative accuracy required in the computed values of the solution.
**thresh** : float
The threshold value for use in the evaluation of the estimated relative errors. For two successive meshes the following condition must hold at each point of the coarser mesh
.. math::
\frac{{\left\lvert Y_1-Y_2\right\rvert }}{{\mathrm{max}\left(\left\lvert Y_1\right\rvert, \left\lvert Y_2\right\rvert, \left\lvert \mathrm{thresh}\right\rvert \right)}}\leq \mathrm{tol}\text{,}
where :math:`Y_1` is the computed solution on the coarser mesh and :math:`Y_2` is the computed solution at the corresponding point in the finer mesh. If this condition is not satisfied then the step size is halved and the solution is recomputed.
Note: :math:`\mathrm{thresh}` can be used to effect a relative, absolute or mixed error test. If :math:`\mathrm{thresh} = 0.0` then pure relative error is measured and, if the computed solution is small and :math:`\mathrm{thresh} = 1.0`, absolute error is measured.
**lwk** : int
The dimension of the array :math:`\mathrm{work}`.
**data** : arbitrary, optional
User-communication data for callback functions.
**Returns**
**yn** : float, ndarray, shape :math:`\left(\mathrm{nmesh}\right)`
:math:`\mathrm{yn}[\textit{i}-1]` contains the most recent approximation of the true solution :math:`y\left(t\right)` at the specified point :math:`t = \mathrm{alim}+\textit{i}\times \textit{H}`, for :math:`\textit{i} = 1,2,\ldots,\mathrm{nmesh}`, where :math:`\textit{H} = \left(\mathrm{tlim}-\mathrm{alim}\right)/\mathrm{nmesh}`.
**errest** : float, ndarray, shape :math:`\left(\mathrm{nmesh}\right)`
:math:`\mathrm{errest}[\textit{i}-1]` contains the most recent approximation of the relative error in the computed solution at the point :math:`t = \mathrm{alim}+\textit{i}\times \textit{H}`, for :math:`\textit{i} = 1,2,\ldots,\mathrm{nmesh}`, where :math:`\textit{H} = \left(\mathrm{tlim}-\mathrm{alim}\right)/\mathrm{nmesh}`.
**work** : float, ndarray, shape :math:`\left(\mathrm{lwk}\right)`
If :math:`\mathrm{errno}` = 5 or 6, :math:`\mathrm{work}[0]` contains the size of :math:`\mathrm{lwk}` required for the algorithm to proceed further.
.. _d05ba-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\mathrm{method} = \texttt{'A'}` and :math:`\mathrm{iorder} = 2`.
Constraint: if :math:`\mathrm{method} = \texttt{'A'}`, :math:`3\leq \mathrm{iorder}\leq 6`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{tol} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\sqrt{\text{machine precision}}\leq \mathrm{tol}\leq 1.0`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{alim} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{tlim} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{tlim} > \mathrm{alim}`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{alim} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{alim}\geq 0.0`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{method} = \texttt{'B'}` and :math:`\mathrm{iorder} = 6`.
Constraint: if :math:`\mathrm{method} = \texttt{'B'}`, :math:`2\leq \mathrm{iorder}\leq 5`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{iorder} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`2\leq \mathrm{iorder}\leq 6`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{method} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{method} = \texttt{'A'}` or :math:`\texttt{'B'}`.
(`errno` :math:`2`)
On entry, :math:`\mathrm{method} = \texttt{'A'}`, :math:`\mathrm{iorder} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{nmesh} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: if :math:`\mathrm{method} = \texttt{'A'}`, :math:`\mathrm{nmesh}\geq \mathrm{iorder}-1`.
(`errno` :math:`2`)
On entry, :math:`\mathrm{method} = \texttt{'B'}`, :math:`\mathrm{iorder} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{nmesh} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: if :math:`\mathrm{method} = \texttt{'B'}`, :math:`\mathrm{nmesh}\geq \mathrm{iorder}`.
(`errno` :math:`3`)
On entry, :math:`\mathrm{lwk} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{lwk}\geq 10\times \mathrm{nmesh}+6`; that is, :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`4`)
The solution is not converging. See `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05baf.html#fcomments>`__.
**Warns**
**NagAlgorithmicWarning**
(`errno` :math:`5`)
The workspace which has been supplied is too small for the required accuracy. The number of extrapolations, so far, is :math:`\langle\mathit{\boldsymbol{value}}\rangle`. If you require one more extrapolation extend the size of workspace to: :math:`\mathrm{lwk} = \langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`6`)
The workspace which has been supplied is too small for the required accuracy. The number of extrapolations, so far, is :math:`\langle\mathit{\boldsymbol{value}}\rangle`. If you require one more extrapolation extend the size of workspace to: :math:`\mathrm{lwk} = \langle\mathit{\boldsymbol{value}}\rangle`.
.. _d05ba-py2-py-notes:
**Notes**
``volterra2`` computes the numerical solution of the nonlinear convolution Volterra integral equation of the second kind
.. math::
y\left(t\right) = f\left(t\right)+\int_a^tk\left(t-s\right)g\left(s, {y\left(s\right)}\right){ds}\text{, }\quad a\leq t\leq T\text{.}
It is assumed that the functions involved in `(1) <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05baf.html#eqn1>`__ are sufficiently smooth.
The function uses a reducible linear multi-step formula selected by you to generate a family of quadrature rules.
The reducible formulae available in ``volterra2`` are the Adams--Moulton formulae of orders :math:`3` to :math:`6`, and the backward differentiation formulae (BDF) of orders :math:`2` to :math:`5`.
For more information about the behaviour and the construction of these rules we refer to Lubich (1983) and Wolkenfelt (1982).
The algorithm is based on computing the solution in a step-by-step fashion on a mesh of equispaced points.
The initial step size which is given by :math:`\left(T-a\right)/N`, :math:`N` being the number of points at which the solution is sought, is halved and another approximation to the solution is computed.
This extrapolation procedure is repeated until successive approximations satisfy a user-specified error requirement.
The above methods require some starting values.
For the Adams' formula of order greater than :math:`3` and the BDF of order greater than :math:`2` we employ an explicit Dormand--Prince--Shampine Runge--Kutta method (see Shampine (1986)).
The above scheme avoids the calculation of the kernel, :math:`k\left(t\right)`, on the negative real line.
.. _d05ba-py2-py-references:
**References**
Lubich, Ch, 1983, `On the stability of linear multi-step methods for Volterra convolution equations`, IMA J. Numer. Anal. (3), 439--465
Shampine, L F, 1986, `Some practical Runge--Kutta formulas`, Math. Comput. (46(173)), 135--150
Wolkenfelt, P H M, 1982, `The construction of reducible quadrature rules for Volterra integral and integro-differential equations`, IMA J. Numer. Anal. (2), 131--152
"""
raise NotImplementedError
[docs]def abel2_weak(ck, cf, cg, initwt, tlim, nmesh, comm, iorder=4, tolnl=sqrt(machine.precision()), data=None):
r"""
``abel2_weak`` computes the solution of a weakly singular nonlinear convolution Volterra--Abel integral equation of the second kind using a fractional Backward Differentiation Formulae (BDF) method.
.. _d05bd-py2-py-doc:
For full information please refer to the NAG Library document for d05bd
https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05bdf.html
.. _d05bd-py2-py-parameters:
**Parameters**
**ck** : callable retval = ck(t, data=None)
:math:`\mathrm{ck}` must evaluate the kernel :math:`k\left(t\right)` of the integral equation `(1) <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05bdf.html#eqn1>`__.
**Parameters**
**t** : float
:math:`t`, the value of the independent variable.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**Returns**
**retval** : float
The value of :math:`k\left(t\right)` evaluated at :math:`\mathrm{t}`.
**cf** : callable retval = cf(t, data=None)
:math:`\mathrm{cf}` must evaluate the function :math:`f\left(t\right)` in `(1) <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05bdf.html#eqn1>`__.
**Parameters**
**t** : float
:math:`t`, the value of the independent variable.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**Returns**
**retval** : float
The value of :math:`f\left(t\right)` evaluated at :math:`\mathrm{t}`.
**cg** : callable retval = cg(s, y, data=None)
:math:`\mathrm{cg}` must evaluate the function :math:`g\left(s, {y\left(s\right)}\right)` in `(1) <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05bdf.html#eqn1>`__.
**Parameters**
**s** : float
:math:`s`, the value of the independent variable.
**y** : float
The value of the solution :math:`y` at the point :math:`\mathrm{s}`.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**Returns**
**retval** : float
The value of :math:`g\left(s, {y\left(s\right)}\right)` evaluated at :math:`\mathrm{s}` and :math:`\mathrm{y}`.
**initwt** : str, length 1
If the fractional weights required by the method need to be calculated by the function then set :math:`\mathrm{initwt} = \texttt{'I'}` (**I**\ nitial call).
If :math:`\mathrm{initwt} = \texttt{'S'}` (**S**\ ubsequent call), the function assumes the fractional weights have been computed on a previous call and are stored in :math:`\mathrm{comm}`\ ['work'].
**tlim** : float
The final point of the integration interval, :math:`\textit{T}`.
**nmesh** : int
:math:`N`, the number of equispaced points at which the solution is sought.
**comm** : dict, communication object, modified in place
Communication structure.
`On initial entry`: need not be set.
**iorder** : int, optional
:math:`p`, the order of the BDF method to be used.
**tolnl** : float, optional
The accuracy required for the computation of the starting value and the solution of the nonlinear equation at each step of the computation (see `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05bdf.html#fcomments>`__).
**data** : arbitrary, optional
User-communication data for callback functions.
**Returns**
**yn** : float, ndarray, shape :math:`\left(\mathrm{nmesh}\right)`
:math:`\mathrm{yn}[\textit{i}-1]` contains the approximate value of the true solution :math:`y\left(t\right)` at the point :math:`t = \left(\textit{i}-1\right)\times h`, for :math:`\textit{i} = 1,2,\ldots,\mathrm{nmesh}`, where :math:`h = \mathrm{tlim}/\left(\mathrm{nmesh}-1\right)`.
.. _d05bd-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\mathrm{tolnl} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{tolnl} > 10\times \text{machine precision}`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{nmesh} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{iorder} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{nmesh} = 2^m+2\times \mathrm{iorder}-1`, for some :math:`m`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{nmesh} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{iorder} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{nmesh}\geq 2\times \mathrm{iorder}+1`.
(`errno` :math:`1`)
On entry, :math:`\textit{lwk} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{lwk}\geq \left(2\times \mathrm{iorder}+6\right)\times \mathrm{nmesh}+8\times \mathrm{iorder}^2-16\times \mathrm{iorder}+1`; that is, :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{tlim} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraints: :math:`\mathrm{tlim} > 10\times \text{machine precision}`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{initwt} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraints: :math:`\mathrm{initwt} = \texttt{'I'}` or :math:`\texttt{'S'}`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{iorder} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`4\leq \mathrm{iorder}\leq 6`.
(`errno` :math:`2`)
An error occurred when trying to compute the starting values.
(`errno` :math:`3`)
An error occurred when trying to compute the solution at a specific step.
.. _d05bd-py2-py-notes:
**Notes**
``abel2_weak`` computes the numerical solution of the weakly singular convolution Volterra--Abel integral equation of the second kind
.. math::
y\left(t\right) = f\left(t\right)+\frac{1}{\sqrt{\pi }}\int_0^t\frac{{k\left(t-s\right)}}{{\sqrt{t-s}}}g\left(s, {y\left(s\right)}\right){ds}\text{, }\quad 0\leq t\leq T\text{.}
Note the constant :math:`\frac{1}{\sqrt{\pi }}` in `(1) <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05bdf.html#eqn1>`__.
It is assumed that the functions involved in `(1) <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05bdf.html#eqn1>`__ are sufficiently smooth.
The function uses a fractional BDF linear multi-step method to generate a family of quadrature rules (see :meth:`abel_weak_weights`).
The BDF methods available in ``abel2_weak`` are of orders :math:`4`, :math:`5` and :math:`6` (:math:`\text{} = p` say).
For a description of the theoretical and practical background to these methods we refer to Lubich (1985) and to Baker and Derakhshan (1987) and Hairer `et al.` (1988) respectively.
The algorithm is based on computing the solution :math:`y\left(t\right)` in a step-by-step fashion on a mesh of equispaced points.
The size of the mesh is given by :math:`T/\left(N-1\right)`, :math:`N` being the number of points at which the solution is sought.
These methods require :math:`2p-1` (including :math:`y\left(0\right)`) starting values which are evaluated internally.
The computation of the lag term arising from the discretization of `(1) <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05bdf.html#eqn1>`__ is performed by fast Fourier transform (FFT) techniques when :math:`N > 32+2p-1`, and directly otherwise.
The function does not provide an error estimate and you are advised to check the behaviour of the solution with a different value of :math:`N`.
An option is provided which avoids the re-evaluation of the fractional weights when ``abel2_weak`` is to be called several times (with the same value of :math:`N`) within the same program unit with different functions.
.. _d05bd-py2-py-references:
**References**
Baker, C T H and Derakhshan, M S, 1987, `FFT techniques in the numerical solution of convolution equations`, J. Comput. Appl. Math. (20), 5--24
Hairer, E, Lubich, Ch and Schlichte, M, 1988, `Fast numerical solution of weakly singular Volterra integral equations`, J. Comput. Appl. Math. (23), 87--98
Lubich, Ch, 1985, `Fractional linear multistep methods for Abel--Volterra integral equations of the second kind`, Math. Comput. (45), 463--469
"""
raise NotImplementedError
[docs]def abel1_weak(ck, cf, cg, initwt, tlim, yn, comm, iorder=4, tolnl=sqrt(machine.precision()), data=None):
r"""
``abel1_weak`` computes the solution of a weakly singular nonlinear convolution Volterra--Abel integral equation of the first kind using a fractional Backward Differentiation Formulae (BDF) method.
.. _d05be-py2-py-doc:
For full information please refer to the NAG Library document for d05be
https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05bef.html
.. _d05be-py2-py-parameters:
**Parameters**
**ck** : callable retval = ck(t, data=None)
:math:`\mathrm{ck}` must evaluate the kernel :math:`k\left(t\right)` of the integral equation `(1) <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05bef.html#eqn1>`__.
**Parameters**
**t** : float
:math:`t`, the value of the independent variable.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**Returns**
**retval** : float
The value of :math:`k\left(t\right)` evaluated at :math:`\mathrm{t}`.
**cf** : callable retval = cf(t, data=None)
:math:`\mathrm{cf}` must evaluate the function :math:`f\left(t\right)` in `(1) <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05bef.html#eqn1>`__.
**Parameters**
**t** : float
:math:`t`, the value of the independent variable.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**Returns**
**retval** : float
The value of :math:`f\left(t\right)` evaluated at :math:`\mathrm{t}`.
**cg** : callable retval = cg(s, y, data=None)
:math:`\mathrm{cg}` must evaluate the function :math:`g\left(s, {y\left(s\right)}\right)` in `(1) <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05bef.html#eqn1>`__.
**Parameters**
**s** : float
:math:`s`, the value of the independent variable.
**y** : float
The value of the solution :math:`y` at the point :math:`\mathrm{s}`.
**data** : arbitrary, optional, modifiable in place
User-communication data for callback functions.
**Returns**
**retval** : float
The value of :math:`g\left(s, {y\left(s\right)}\right)` evaluated at :math:`\mathrm{s}` and :math:`\mathrm{y}`.
**initwt** : str, length 1
If the fractional weights required by the method need to be calculated by the function then set :math:`\mathrm{initwt} = \texttt{'I'}` (**I**\ nitial call).
If :math:`\mathrm{initwt} = \texttt{'S'}` (**S**\ ubsequent call), the function assumes the fractional weights have been computed by a previous call and are stored in :math:`\mathrm{comm}`\ ['work'].
**tlim** : float
The final point of the integration interval, :math:`\textit{T}`.
**yn** : float, array-like, shape :math:`\left(\textit{nmesh}\right)`
:math:`\mathrm{yn}[0]` must contain the value of :math:`y\left(t\right)` at :math:`t = 0` (see `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05bef.html#fcomments>`__).
**comm** : dict, communication object, modified in place
Communication structure.
`On initial entry`: need not be set.
**iorder** : int, optional
:math:`p`, the order of the BDF method to be used.
**tolnl** : float, optional
The accuracy required for the computation of the starting value and the solution of the nonlinear equation at each step of the computation (see `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05bef.html#fcomments>`__).
**data** : arbitrary, optional
User-communication data for callback functions.
**Returns**
**yn** : float, ndarray, shape :math:`\left(\textit{nmesh}\right)`
:math:`\mathrm{yn}[\textit{i}-1]` contains the approximate value of the true solution :math:`y\left(t\right)` at the point :math:`t = \left(\textit{i}-1\right)\times h`, for :math:`\textit{i} = 1,2,\ldots,\textit{nmesh}`, where :math:`h = \mathrm{tlim}/\left(\textit{nmesh}-1\right)`.
.. _d05be-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\mathrm{tolnl} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{tolnl} > 10\times \text{machine precision}`.
(`errno` :math:`1`)
On entry, :math:`\textit{nmesh} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{iorder} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nmesh} = 2^m+2\times \mathrm{iorder}-1`, for some :math:`m`.
(`errno` :math:`1`)
On entry, :math:`\textit{nmesh} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{iorder} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{nmesh}\geq 2\times \mathrm{iorder}+1`.
(`errno` :math:`1`)
On entry, :math:`\textit{lwk} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{lwk}\geq \left(2\times \mathrm{iorder}+6\right)\times \textit{nmesh}+8\times \mathrm{iorder}^2-16\times \mathrm{iorder}+1`; that is, :math:`\langle\mathit{\boldsymbol{value}}\rangle`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{tlim} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{tlim} > 10\times \text{machine precision}`
(`errno` :math:`1`)
On entry, :math:`\mathrm{initwt} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{initwt} = \texttt{'I'}` or :math:`\texttt{'S'}`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{iorder} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`4\leq \mathrm{iorder}\leq 6`.
(`errno` :math:`2`)
An error occurred when trying to compute the starting values.
(`errno` :math:`3`)
An error occurred when trying to compute the solution at a specific step.
.. _d05be-py2-py-notes:
**Notes**
``abel1_weak`` computes the numerical solution of the weakly singular convolution Volterra--Abel integral equation of the first kind
.. math::
f\left(t\right)+\frac{1}{\sqrt{\pi }}\int_0^t\frac{{k\left(t-s\right)}}{{\sqrt{t-s}}}g\left(s, {y\left(s\right)}\right){ds} = 0\text{, }\quad 0\leq t\leq T\text{.}
Note the constant :math:`\frac{1}{\sqrt{\pi }}` in `(1) <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05bef.html#eqn1>`__.
It is assumed that the functions involved in `(1) <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05bef.html#eqn1>`__ are sufficiently smooth and if
.. math::
f\left(t\right) = t^{\beta }w\left(t\right)\quad \text{ with }\quad \beta > -\frac{1}{2}\text{ and }w\left(t\right)\text{ smooth,}
then the solution :math:`y\left(t\right)` is unique and has the form :math:`y\left(t\right) = t^{{\beta -1/2}}z\left(t\right)`, (see Lubich (1987)).
It is evident from `(1) <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05bef.html#eqn1>`__ that :math:`f\left(0\right) = 0`.
You are required to provide the value of :math:`y\left(t\right)` at :math:`t = 0`.
If :math:`y\left(0\right)` is unknown, `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05bef.html#fcomments>`__ gives a description of how an approximate value can be obtained.
The function uses a fractional BDF linear multi-step method selected by you to generate a family of quadrature rules (see :meth:`abel_weak_weights`).
The BDF methods available in ``abel1_weak`` are of orders :math:`4`, :math:`5` and :math:`6` (:math:`\text{} = p` say).
For a description of the theoretical and practical background related to these methods we refer to Lubich (1987) and to Baker and Derakhshan (1987) and Hairer `et al.` (1988) respectively.
The algorithm is based on computing the solution :math:`y\left(t\right)` in a step-by-step fashion on a mesh of equispaced points.
The size of the mesh is given by :math:`T/\left(N-1\right)`, :math:`N` being the number of points at which the solution is sought.
These methods require :math:`2p-2` starting values which are evaluated internally.
The computation of the lag term arising from the discretization of `(1) <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05bef.html#eqn1>`__ is performed by fast Fourier transform (FFT) techniques when :math:`N > 32+2p-1`, and directly otherwise.
The function does not provide an error estimate and you are advised to check the behaviour of the solution with a different value of :math:`N`.
An option is provided which avoids the re-evaluation of the fractional weights when ``abel1_weak`` is to be called several times (with the same value of :math:`N`) within the same program with different functions.
.. _d05be-py2-py-references:
**References**
Baker, C T H and Derakhshan, M S, 1987, `FFT techniques in the numerical solution of convolution equations`, J. Comput. Appl. Math. (20), 5--24
Gorenflo, R and Pfeiffer, A, 1991, `On analysis and discretization of nonlinear Abel integral equations of first kind`, Acta Math. Vietnam (16), 211--262
Hairer, E, Lubich, Ch and Schlichte, M, 1988, `Fast numerical solution of weakly singular Volterra integral equations`, J. Comput. Appl. Math. (23), 87--98
Lubich, Ch, 1987, `Fractional linear multistep methods for Abel--Volterra integral equations of the first kind`, IMA J. Numer. Anal (7), 97--106
"""
raise NotImplementedError
[docs]def volterra_weights(method, iorder, nomg, nwt):
r"""
``volterra_weights`` computes the quadrature weights associated with the Adams' methods of orders three to six and the Backward Differentiation Formulae (BDF) methods of orders two to five.
These rules, which are referred to as reducible quadrature rules, can then be used in the solution of Volterra integral and integro-differential equations.
.. _d05bw-py2-py-doc:
For full information please refer to the NAG Library document for d05bw
https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05bwf.html
.. _d05bw-py2-py-parameters:
**Parameters**
**method** : str, length 1
The type of method to be used.
:math:`\mathrm{method} = \texttt{'A'}`
For Adams' type formulae.
:math:`\mathrm{method} = \texttt{'B'}`
For Backward Differentiation Formulae.
**iorder** : int
The order of the method to be used. The number of starting weights, :math:`\textit{p}` is determined by :math:`\mathrm{method}` and :math:`\mathrm{iorder}`.
If :math:`\mathrm{method} = \texttt{'A'}`, :math:`\textit{p} = \mathrm{iorder}-1`.
If :math:`\mathrm{method} = \texttt{'B'}`, :math:`\textit{p} = \mathrm{iorder}`.
**nomg** : int
The number of convolution weights, :math:`\textit{n}_w`.
**nwt** : int
:math:`\textit{p}`, the number of columns in the starting weights.
**Returns**
**omega** : float, ndarray, shape :math:`\left(\mathrm{nomg}\right)`
Contains the first :math:`\mathrm{nomg}` convolution weights.
**lensw** : int
The number of rows in the weights :math:`W_{{i,j}}`.
**sw** : float, ndarray, shape :math:`\left(\textit{n}, \mathrm{nwt}\right)`
:math:`\mathrm{sw}[\textit{i}-1,\textit{j}]` contains the weights :math:`W_{{\textit{i},\textit{j}}}`, for :math:`\textit{j} = 0,1,\ldots,\mathrm{nwt}-1`, for :math:`\textit{i} = 1,2,\ldots,\mathrm{lensw}`, where :math:`\textit{n}` is as defined in :ref:`Notes <d05bw-py2-py-notes>`.
.. _d05bw-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\mathrm{method} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{method} = \texttt{'A'}` or :math:`\texttt{'B'}`.
(`errno` :math:`2`)
On entry, :math:`\mathrm{nomg} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{nomg}\geq 1`.
(`errno` :math:`2`)
On entry, :math:`\mathrm{iorder} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`2\leq \mathrm{iorder}\leq 6`.
(`errno` :math:`3`)
On entry, :math:`\mathrm{method} = \texttt{'B'}` and :math:`\mathrm{iorder} = 6`.
Constraint: if :math:`\mathrm{method} = \texttt{'B'}`, :math:`2\leq \mathrm{iorder}\leq 5`.
(`errno` :math:`3`)
On entry, :math:`\mathrm{method} = \texttt{'A'}` and :math:`\mathrm{iorder} = 2`.
Constraint: if :math:`\mathrm{method} = \texttt{'A'}`, :math:`3\leq \mathrm{iorder}\leq 6`.
(`errno` :math:`4`)
On entry, :math:`\mathrm{method} = \texttt{'B'}`, :math:`\mathrm{iorder} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{nwt} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: if :math:`\mathrm{method} = \texttt{'B'}`, :math:`\mathrm{nwt} = \mathrm{iorder}`.
(`errno` :math:`4`)
On entry, :math:`\mathrm{method} = \texttt{'A'}`, :math:`\mathrm{iorder} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{nwt} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: if :math:`\mathrm{method} = \texttt{'A'}`, :math:`\mathrm{nwt} = \mathrm{iorder}-1`.
(`errno` :math:`5`)
On entry, :math:`\mathrm{method} = \texttt{'B'}`, :math:`\mathrm{iorder} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{nomg} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\textit{ldsw} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: if :math:`\mathrm{method} = \texttt{'B'}`, :math:`\textit{ldsw}\geq \mathrm{nomg}+\mathrm{iorder}-1`.
(`errno` :math:`5`)
On entry, :math:`\mathrm{method} = \texttt{'A'}`, :math:`\mathrm{iorder} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{nomg} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\textit{ldsw} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: if :math:`\mathrm{method} = \texttt{'A'}`, :math:`\textit{ldsw}\geq \mathrm{nomg}+\mathrm{iorder}-2`.
.. _d05bw-py2-py-notes:
**Notes**
``volterra_weights`` computes the weights :math:`W_{{i,j}}` and :math:`\omega_i` for a family of quadrature rules related to the Adams' methods of orders three to six and the BDF methods of orders two to five, for approximating the integral:
.. math::
\int_0^t\phi \left(s\right){ds}\simeq h\sum_{{j = 0}}^{{\textit{p}-1}}W_{{i,j}}\phi \left(j\times h\right)+h\sum_{{j = \textit{p}}}^i\omega_{{i-j}}\phi \left(j\times h\right)\text{, }\quad 0\leq t\leq T\text{,}
with :math:`t = \textit{i}\times h`, for :math:`\textit{i} = 0,1,\ldots,\textit{n}`, for some given constant :math:`h`.
In `(1) <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05bwf.html#eqn1>`__, :math:`h` is a uniform mesh, :math:`\textit{p}` is related to the order of the method being used and :math:`W_{{i,j}}`, :math:`\omega_i` are the starting and the convolution weights respectively.
The mesh size :math:`h` is determined as :math:`h = \frac{T}{\textit{n}}`, where :math:`\textit{n} = \textit{n}_w+\textit{p}-1` and :math:`\textit{n}_w` is the chosen number of convolution weights :math:`w_j`, for :math:`\textit{j} = 1,2,\ldots,{\textit{n}_w-1}`.
A description of how these weights can be used in the solution of a Volterra integral equation of the second kind is given in `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05bwf.html#fcomments>`__.
For a general discussion of these methods, see Wolkenfelt (1982) for more details.
.. _d05bw-py2-py-references:
**References**
Lambert, J D, 1973, `Computational Methods in Ordinary Differential Equations`, John Wiley
Wolkenfelt, P H M, 1982, `The construction of reducible quadrature rules for Volterra integral and integro-differential equations`, IMA J. Numer. Anal. (2), 131--152
"""
raise NotImplementedError
[docs]def abel_weak_weights(iorder, iq, lenfw):
r"""
``abel_weak_weights`` computes the fractional quadrature weights associated with the Backward Differentiation Formulae (BDF) of orders :math:`4`, :math:`5` and :math:`6`.
These weights can then be used in the solution of weakly singular equations of Abel type.
.. _d05by-py2-py-doc:
For full information please refer to the NAG Library document for d05by
https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05byf.html
.. _d05by-py2-py-parameters:
**Parameters**
**iorder** : int
:math:`p`, the order of the BDF method to be used.
**iq** : int
Determines the number of weights to be computed. By setting :math:`\mathrm{iq}` to a value, :math:`2^{{\mathrm{iq}+1}}` fractional convolution weights are computed.
**lenfw** : int
The dimension of the array :math:`\mathrm{wt}`.
**Returns**
**wt** : float, ndarray, shape :math:`\left(\mathrm{lenfw}\right)`
The first :math:`2^{{\mathrm{iq}+1}}` elements of :math:`\mathrm{wt}` contains the fractional convolution weights :math:`\omega_i`, for :math:`\textit{i} = 0,1,\ldots,2^{{\mathrm{iq}+1}}-1`. The remainder of the array is used as workspace.
**sw** : float, ndarray, shape :math:`\left(2^{{\mathrm{iq}+1}}+2\times \mathrm{iorder}-1, 2\times \mathrm{iorder}-1\right)`
:math:`\mathrm{sw}[\textit{i}-1,\textit{j}]` contains the fractional starting weights :math:`W_{{\textit{i}-1,\textit{j}}}`, for :math:`\textit{j} = 0,1,\ldots,2\times \mathrm{iorder}-2`, for :math:`\textit{i} = 1,2,\ldots,\textit{N}`, where :math:`\textit{N} = \left(2^{{\mathrm{iq}+1}}+2\times \mathrm{iorder}-1\right)`.
.. _d05by-py2-py-errors:
**Raises**
**NagValueError**
(`errno` :math:`1`)
On entry, :math:`\mathrm{lenfw} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{iq} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{lenfw}\geq 2^{{\mathrm{iq}+2}}`.
(`errno` :math:`1`)
On entry, :math:`\textit{ldsw} = \langle\mathit{\boldsymbol{value}}\rangle`, :math:`\mathrm{iq} = \langle\mathit{\boldsymbol{value}}\rangle` and :math:`\mathrm{iorder} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\textit{ldsw}\geq 2^{{\mathrm{iq}+1}}+2\times \mathrm{iorder}-1`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{iq} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`\mathrm{iq}\geq 0`.
(`errno` :math:`1`)
On entry, :math:`\mathrm{iorder} = \langle\mathit{\boldsymbol{value}}\rangle`.
Constraint: :math:`4\leq \mathrm{iorder}\leq 6`.
.. _d05by-py2-py-notes:
**Notes**
``abel_weak_weights`` computes the weights :math:`W_{{i,j}}` and :math:`\omega_i` for a family of quadrature rules related to a BDF method for approximating the integral:
.. math::
\frac{1}{\sqrt{\pi }}\int_0^t\frac{{\phi \left(s\right)}}{{\sqrt{t-s}}}{ds}\simeq \sqrt{h}\sum_{{j = 0}}^{{2p-2}}W_{{i,j}}\phi \left(j\times h\right)+\sqrt{h}\sum_{{j = 2p-1}}^i\omega_{{i-j}}\phi \left(j\times h\right)\text{, }\quad 0\leq t\leq T\text{,}
with :math:`t = i\times h\left(i\geq 0\right)`, for some given :math:`h`.
In `(1) <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05byf.html#eqn1>`__, :math:`p` is the order of the BDF method used and :math:`W_{{i,j}}`, :math:`\omega_i` are the fractional starting and the fractional convolution weights respectively.
The algorithm for the generation of :math:`\omega_i` is based on Newton's iteration.
Fast Fourier transform (FFT) techniques are used for computing these weights and subsequently :math:`W_{{i,j}}` (see Baker and Derakhshan (1987) and Henrici (1979) for practical details and Lubich (1986) for theoretical details).
Some special functions can be represented as the fractional integrals of simpler functions and fractional quadratures can be employed for their computation (see Lubich (1986)).
A description of how these weights can be used in the solution of weakly singular equations of Abel type is given in `Further Comments <https://support.nag.com/numeric/nl/nagdoc_30.2/flhtml/d05/d05byf.html#fcomments>`__.
.. _d05by-py2-py-references:
**References**
Baker, C T H and Derakhshan, M S, 1987, `Computational approximations to some power series`, Approximation Theory, (eds L Collatz, G Meinardus and G Nürnberger) (81), 11--20
Henrici, P, 1979, `Fast Fourier methods in computational complex analysis`, SIAM Rev. (21), 481--529
Lubich, Ch, 1986, `Discretized fractional calculus`, SIAM J. Math. Anal. (17), 704--719
"""
raise NotImplementedError
