NAG Library Function Document
nag_inteq_volterra2 (d05bac)
1
Purpose
nag_inteq_volterra2 (d05bac) computes the solution of a nonlinear convolution Volterra integral equation of the second kind using a reducible linear multistep method.
2
Specification
#include <nag.h> 
#include <nagd05.h> 
void 
nag_inteq_volterra2 (
double 
(*ck)(double t,
Nag_Comm *comm),


double 
(*cg)(double s,
double y,
Nag_Comm *comm),


double 
(*cf)(double t,
Nag_Comm *comm),


Nag_ODEMethod method,
Integer iorder,
double alim,
double tlim,
double tol,
Integer nmesh,
double thresh,
double work[],
Integer lwk,
double yn[],
double errest[],
Nag_Comm *comm,
NagError *fail) 

3
Description
nag_inteq_volterra2 (d05bac) computes the numerical solution of the nonlinear convolution Volterra integral equation of the second kind
It is assumed that the functions involved in
(1) are sufficiently smooth. The function uses a reducible linear multistep formula selected by you to generate a family of quadrature rules. The reducible formulae available in
nag_inteq_volterra2 (d05bac) are the Adams–Moulton formulae of orders
$3$ to
$6$, and the backward differentiation formulae (BDF) of orders
$2$ to
$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 stepbystep fashion on a mesh of equispaced points. The initial step size which is given by $\left(Ta\right)/N$, $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 userspecified error requirement.
The above methods require some starting values. For the Adams' formula of order greater than
$3$ and the BDF of order greater than
$2$ we employ an explicit Dormand–Prince–Shampine Runge–Kutta method (see
Shampine (1986)). The above scheme avoids the calculation of the kernel,
$k\left(t\right)$, on the negative real line.
4
References
Lubich Ch (1983) On the stability of linear multistep 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 integrodifferential equations IMA J. Numer. Anal. 2 131–152
5
Arguments
 1:
$\mathbf{ck}$ – function, supplied by the userExternal Function

ck must evaluate the kernel
$k\left(t\right)$ of the integral equation
(1).
The specification of
ck is:
double 
ck (double t,
Nag_Comm *comm)


 1:
$\mathbf{t}$ – doubleInput

On entry: $t$, the value of the independent variable.
 2:
$\mathbf{comm}$ – Nag_Comm *
Pointer to structure of type Nag_Comm; the following members are relevant to
ck.
 user – double *
 iuser – Integer *
 p – Pointer
The type Pointer will be
void *. Before calling
nag_inteq_volterra2 (d05bac) you may allocate memory and initialize these pointers with various quantities for use by
ck when called from
nag_inteq_volterra2 (d05bac) (see
Section 3.3.1.1 in How to Use the NAG Library and its Documentation).
Note: ck should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
nag_inteq_volterra2 (d05bac). If your code inadvertently
does return any NaNs or infinities,
nag_inteq_volterra2 (d05bac) is likely to produce unexpected results.
 2:
$\mathbf{cg}$ – function, supplied by the userExternal Function

cg must evaluate the function
$g\left(s,y\left(s\right)\right)$ in
(1).
The specification of
cg is:
double 
cg (double s,
double y,
Nag_Comm *comm)


 1:
$\mathbf{s}$ – doubleInput

On entry: $s$, the value of the independent variable.
 2:
$\mathbf{y}$ – doubleInput

On entry: the value of the solution
$y$ at the point
s.
 3:
$\mathbf{comm}$ – Nag_Comm *
Pointer to structure of type Nag_Comm; the following members are relevant to
cg.
 user – double *
 iuser – Integer *
 p – Pointer
The type Pointer will be
void *. Before calling
nag_inteq_volterra2 (d05bac) you may allocate memory and initialize these pointers with various quantities for use by
cg when called from
nag_inteq_volterra2 (d05bac) (see
Section 3.3.1.1 in How to Use the NAG Library and its Documentation).
Note: cg should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
nag_inteq_volterra2 (d05bac). If your code inadvertently
does return any NaNs or infinities,
nag_inteq_volterra2 (d05bac) is likely to produce unexpected results.
 3:
$\mathbf{cf}$ – function, supplied by the userExternal Function

cf must evaluate the function
$f\left(t\right)$ in
(1).
The specification of
cf is:
double 
cf (double t,
Nag_Comm *comm)


 1:
$\mathbf{t}$ – doubleInput

On entry: $t$, the value of the independent variable.
 2:
$\mathbf{comm}$ – Nag_Comm *
Pointer to structure of type Nag_Comm; the following members are relevant to
cf.
 user – double *
 iuser – Integer *
 p – Pointer
The type Pointer will be
void *. Before calling
nag_inteq_volterra2 (d05bac) you may allocate memory and initialize these pointers with various quantities for use by
cf when called from
nag_inteq_volterra2 (d05bac) (see
Section 3.3.1.1 in How to Use the NAG Library and its Documentation).
Note: cf should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
nag_inteq_volterra2 (d05bac). If your code inadvertently
does return any NaNs or infinities,
nag_inteq_volterra2 (d05bac) is likely to produce unexpected results.
 4:
$\mathbf{method}$ – Nag_ODEMethodInput

On entry: the type of method which you wish to employ.
 ${\mathbf{method}}=\mathrm{Nag\_Adams}$
 For Adams' type formulae.
 ${\mathbf{method}}=\mathrm{Nag\_BDF}$
 For backward differentiation formulae.
Constraint:
${\mathbf{method}}=\mathrm{Nag\_Adams}$ or $\mathrm{Nag\_BDF}$.
 5:
$\mathbf{iorder}$ – IntegerInput

On entry: the order of the method to be used.
Constraints:
 if ${\mathbf{method}}=\mathrm{Nag\_Adams}$, $3\le {\mathbf{iorder}}\le 6$;
 if ${\mathbf{method}}=\mathrm{Nag\_BDF}$, $2\le {\mathbf{iorder}}\le 5$.
 6:
$\mathbf{alim}$ – doubleInput

On entry: $a$, the lower limit of the integration interval.
Constraint:
${\mathbf{alim}}\ge 0.0$.
 7:
$\mathbf{tlim}$ – doubleInput

On entry: the final point of the integration interval, $T$.
Constraint:
${\mathbf{tlim}}>{\mathbf{alim}}$.
 8:
$\mathbf{tol}$ – doubleInput

On entry: the relative accuracy required in the computed values of the solution.
Constraint:
$\sqrt{\epsilon}\le {\mathbf{tol}}\le 1.0$, where $\epsilon $ is the machine precision.
 9:
$\mathbf{nmesh}$ – IntegerInput

On entry: the number of equidistant points at which the solution is sought.
Constraints:
 if ${\mathbf{method}}=\mathrm{Nag\_Adams}$, ${\mathbf{nmesh}}\ge {\mathbf{iorder}}1$;
 if ${\mathbf{method}}=\mathrm{Nag\_BDF}$, ${\mathbf{nmesh}}\ge {\mathbf{iorder}}$.
 10:
$\mathbf{thresh}$ – doubleInput

On entry: 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
where
${Y}_{1}$ is the computed solution on the coarser mesh and
${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: thresh can be used to effect a relative, absolute or mixed error test. If
${\mathbf{thresh}}=0.0$ then pure relative error is measured and, if the computed solution is small and
${\mathbf{thresh}}=1.0$, absolute error is measured.
 11:
$\mathbf{work}\left[{\mathbf{lwk}}\right]$ – doubleOutput
 12:
$\mathbf{lwk}$ – IntegerInput

On entry: the dimension of the array
work.
Constraint:
${\mathbf{lwk}}\ge 10\times {\mathbf{nmesh}}+6$.
Note: the above value of
lwk is sufficient for
nag_inteq_volterra2 (d05bac) to perform only one extrapolation on the initial mesh as defined by
nmesh. In general much more workspace is required and in the case when a large step size is supplied (i.e.,
nmesh is small), you must provide a considerably larger workspace.
On exit: if
${\mathbf{fail}}\mathbf{.}\mathbf{code}=$ NW_OUT_OF_WORKSPACE,
${\mathbf{work}}\left[0\right]$ contains the size of
lwk required for the algorithm to proceed further.
 13:
$\mathbf{yn}\left[{\mathbf{nmesh}}\right]$ – doubleOutput

On exit: ${\mathbf{yn}}\left[\mathit{i}1\right]$ contains the most recent approximation of the true solution $y\left(t\right)$ at the specified point $t={\mathbf{alim}}+\mathit{i}\times H$, for $\mathit{i}=1,2,\dots ,{\mathbf{nmesh}}$, where $H=\left({\mathbf{tlim}}{\mathbf{alim}}\right)/{\mathbf{nmesh}}$.
 14:
$\mathbf{errest}\left[{\mathbf{nmesh}}\right]$ – doubleOutput

On exit: ${\mathbf{errest}}\left[\mathit{i}1\right]$ contains the most recent approximation of the relative error in the computed solution at the point $t={\mathbf{alim}}+\mathit{i}\times H$, for $\mathit{i}=1,2,\dots ,{\mathbf{nmesh}}$, where $H=\left({\mathbf{tlim}}{\mathbf{alim}}\right)/{\mathbf{nmesh}}$.
 15:
$\mathbf{comm}$ – Nag_Comm *

The NAG communication argument (see
Section 3.3.1.1 in How to Use the NAG Library and its Documentation).
 16:
$\mathbf{fail}$ – NagError *Input/Output

The NAG error argument (see
Section 3.7 in How to Use the NAG Library and its Documentation).
6
Error Indicators and Warnings
 NE_ALLOC_FAIL

Dynamic memory allocation failed.
See
Section 2.3.1.2 in How to Use the NAG Library and its Documentation for further information.
 NE_BAD_PARAM

On entry, argument $\u2329\mathit{\text{value}}\u232a$ had an illegal value.
 NE_CONVERGENCE

The solution is not converging. See
Section 9.
 NE_ENUM_INT

On entry, ${\mathbf{method}}=\mathrm{Nag\_Adams}$ and ${\mathbf{iorder}}=2$.
Constraint: if ${\mathbf{method}}=\mathrm{Nag\_Adams}$, $3\le {\mathbf{iorder}}\le 6$.
On entry, ${\mathbf{method}}=\mathrm{Nag\_BDF}$ and ${\mathbf{iorder}}=6$.
Constraint: if ${\mathbf{method}}=\mathrm{Nag\_BDF}$, $2\le {\mathbf{iorder}}\le 5$.
 NE_ENUM_INT_2

On entry, ${\mathbf{method}}=\mathrm{Nag\_Adams}$, ${\mathbf{iorder}}=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{nmesh}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: if ${\mathbf{method}}=\mathrm{Nag\_Adams}$, ${\mathbf{nmesh}}\ge {\mathbf{iorder}}1$.
On entry, ${\mathbf{method}}=\mathrm{Nag\_BDF}$, ${\mathbf{iorder}}=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{nmesh}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: if ${\mathbf{method}}=\mathrm{Nag\_BDF}$, ${\mathbf{nmesh}}\ge {\mathbf{iorder}}$.
 NE_INT

On entry, ${\mathbf{iorder}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: $2\le {\mathbf{iorder}}\le 6$.
On entry, ${\mathbf{lwk}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{lwk}}\ge 10\times {\mathbf{nmesh}}+6$; that is, $\u2329\mathit{\text{value}}\u232a$.
 NE_INTERNAL_ERROR

An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact
NAG for assistance.
See
Section 2.7.6 in How to Use the NAG Library and its Documentation for further information.
 NE_NO_LICENCE

Your licence key may have expired or may not have been installed correctly.
See
Section 2.7.5 in How to Use the NAG Library and its Documentation for further information.
 NE_REAL

On entry, ${\mathbf{alim}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{alim}}\ge 0.0$.
On entry, ${\mathbf{tol}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: $\sqrt{\mathit{machineprecision}}\le {\mathbf{tol}}\le 1.0$.
 NE_REAL_2

On entry, ${\mathbf{alim}}=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{tlim}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{tlim}}>{\mathbf{alim}}$.
 NW_OUT_OF_WORKSPACE

The workspace which has been supplied is too small for the required accuracy. The number of extrapolations, so far, is $\u2329\mathit{\text{value}}\u232a$. If you require one more extrapolation extend the size of workspace to: ${\mathbf{lwk}}=\u2329\mathit{\text{value}}\u232a$.
7
Accuracy
The accuracy depends on
tol, the theoretical behaviour of the solution of the integral equation, the interval of integration and on the method being used. It can be controlled by varying
tol and
thresh; you are recommended to choose a smaller value for
tol, the larger the value of
iorder.
You are warned not to supply a very small
tol, because the required accuracy may never be achieved. This will usually force an error exit with
${\mathbf{fail}}\mathbf{.}\mathbf{code}=$ NW_OUT_OF_WORKSPACE.
In general, the higher the order of the method, the faster the required accuracy is achieved with less workspace. For nonstiff problems (see
Section 9) you are recommended to use the Adams' method (
${\mathbf{method}}=\mathrm{Nag\_Adams}$) of order greater than
$4$ (
${\mathbf{iorder}}>4$).
8
Parallelism and Performance
nag_inteq_volterra2 (d05bac) makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the
x06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this function. Please also consult the
Users' Note for your implementation for any additional implementationspecific information.
When solving
(1), the solution of a nonlinear equation of the form
is required, where
${\Psi}_{n}$ and
$\alpha $ are constants.
nag_inteq_volterra2 (d05bac) calls
nag_interval_zero_cont_func (c05avc) to find an interval for the zero of this equation followed by
nag_zero_cont_func_brent_rcomm (c05azc) to find its zero.
There is an initial phase of the algorithm where the solution is computed only for the first few points of the mesh. The exact number of these points depends on
iorder and
method. The step size is halved until the accuracy requirements are satisfied on these points and only then the solution on the whole mesh is computed. During this initial phase, if
lwk is too small,
nag_inteq_volterra2 (d05bac) will exit with
${\mathbf{fail}}\mathbf{.}\mathbf{code}=$ NW_OUT_OF_WORKSPACE.
In the case
${\mathbf{fail}}\mathbf{.}\mathbf{code}=$ NE_CONVERGENCE or
NW_OUT_OF_WORKSPACE, you may be dealing with a ‘stiff’ equation; an equation where the Lipschitz constant
$L$ of the function
$g\left(t,y\right)$ in
(1) with respect to its second argument is large, viz,
In this case, if a BDF method (
${\mathbf{method}}=\mathrm{Nag\_BDF}$) has been used, you are recommended to choose a smaller step size by increasing the value of
nmesh, or provide a larger workspace. But, if an Adams' method (
${\mathbf{method}}=\mathrm{Nag\_Adams}$) has been selected, you are recommended to switch to a BDF method instead.
In the case
${\mathbf{fail}}\mathbf{.}\mathbf{code}=$ NW_OUT_OF_WORKSPACE, then if
${\mathbf{fail}}\mathbf{.}\mathbf{errnum}=6$,
the specified accuracy has not been attained but
yn and
errest contain the most recent approximation to the computed solution and the corresponding error estimate. In this case, the error message informs you of the number of extrapolations performed and the size of
lwk required for the algorithm to proceed further. The latter quantity will also be available in
${\mathbf{work}}\left[0\right]$.
10
Example
Consider the following integral equation
with the solution
$y\left(t\right)=\mathrm{ln}\left(t+e\right)$. In this example, the Adams' method of order
$6$ is used to solve this equation with
${\mathbf{tol}}=\text{1.e\u22124}$.
10.1
Program Text
Program Text (d05bace.c)
10.2
Program Data
None.
10.3
Program Results
Program Results (d05bace.r)