NAG CL Interface
d05bac (volterra2)

1 Purpose

d05bac computes the solution of a nonlinear convolution Volterra integral equation of the second kind using a reducible linear multi-step method.

2 Specification

#include <nag.h>
void  d05bac (
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)
The function may be called by the names: d05bac or nag_inteq_volterra2.

3 Description

d05bac computes the numerical solution of the nonlinear convolution Volterra integral equation of the second kind
yt=ft+atkt-sgs,ysds,  atT. (1)
It is assumed that the functions involved in (1) 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 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 step-by-step fashion on a mesh of equispaced points. The initial step size which is given by T-a/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 user-specified 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, kt, on the negative real line.

4 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

5 Arguments

1: ck function, supplied by the user External Function
ck must evaluate the kernel kt of the integral equation (1).
The specification of ck is:
double  ck (double t, Nag_Comm *comm)
1: t double Input
On entry: t, the value of the independent variable.
2: comm Nag_Comm *
Pointer to structure of type Nag_Comm; the following members are relevant to ck.
userdouble *
iuserInteger *
pPointer 
The type Pointer will be void *. Before calling d05bac you may allocate memory and initialize these pointers with various quantities for use by ck when called from d05bac (see Section 3.1.1 in the Introduction to the NAG Library CL Interface).
Note: ck should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d05bac. If your code inadvertently does return any NaNs or infinities, d05bac is likely to produce unexpected results.
2: cg function, supplied by the user External Function
cg must evaluate the function gs,ys in (1).
The specification of cg is:
double  cg (double s, double y, Nag_Comm *comm)
1: s double Input
On entry: s, the value of the independent variable.
2: y double Input
On entry: the value of the solution y at the point s.
3: comm Nag_Comm *
Pointer to structure of type Nag_Comm; the following members are relevant to cg.
userdouble *
iuserInteger *
pPointer 
The type Pointer will be void *. Before calling d05bac you may allocate memory and initialize these pointers with various quantities for use by cg when called from d05bac (see Section 3.1.1 in the Introduction to the NAG Library CL Interface).
Note: cg should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d05bac. If your code inadvertently does return any NaNs or infinities, d05bac is likely to produce unexpected results.
3: cf function, supplied by the user External Function
cf must evaluate the function ft in (1).
The specification of cf is:
double  cf (double t, Nag_Comm *comm)
1: t double Input
On entry: t, the value of the independent variable.
2: comm Nag_Comm *
Pointer to structure of type Nag_Comm; the following members are relevant to cf.
userdouble *
iuserInteger *
pPointer 
The type Pointer will be void *. Before calling d05bac you may allocate memory and initialize these pointers with various quantities for use by cf when called from d05bac (see Section 3.1.1 in the Introduction to the NAG Library CL Interface).
Note: cf should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d05bac. If your code inadvertently does return any NaNs or infinities, d05bac is likely to produce unexpected results.
4: method Nag_ODEMethod Input
On entry: the type of method which you wish to employ.
method=Nag_Adams
For Adams' type formulae.
method=Nag_BDF
For backward differentiation formulae.
Constraint: method=Nag_Adams or Nag_BDF.
5: iorder Integer Input
On entry: the order of the method to be used.
Constraints:
  • if method=Nag_Adams, 3iorder6;
  • if method=Nag_BDF, 2iorder5.
6: alim double Input
On entry: a, the lower limit of the integration interval.
Constraint: alim0.0.
7: tlim double Input
On entry: the final point of the integration interval, T.
Constraint: tlim>alim.
8: tol double Input
On entry: the relative accuracy required in the computed values of the solution.
Constraint: εtol1.0, where ε is the machine precision.
9: nmesh Integer Input
On entry: the number of equidistant points at which the solution is sought.
Constraints:
  • if method=Nag_Adams, nmeshiorder-1;
  • if method=Nag_BDF, nmeshiorder.
10: thresh double Input
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
Y1-Y2 maxY1,Y2,thresh tol,  
where Y1 is the computed solution on the coarser mesh and Y2 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 thresh=0.0 then pure relative error is measured and, if the computed solution is small and thresh=1.0, absolute error is measured.
11: work[lwk] double Output
12: lwk Integer Input
On entry: the dimension of the array work.
Constraint: lwk10×nmesh+6.
Note: the above value of lwk is sufficient for 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 fail.code= NW_OUT_OF_WORKSPACE, work[0] contains the size of lwk required for the algorithm to proceed further.
13: yn[nmesh] double Output
On exit: yn[i-1] contains the most recent approximation of the true solution yt at the specified point t=alim+i×H, for i=1,2,,nmesh, where H=tlim-alim/nmesh.
14: errest[nmesh] double Output
On exit: errest[i-1] contains the most recent approximation of the relative error in the computed solution at the point t=alim+i×H, for i=1,2,,nmesh, where H=tlim-alim/nmesh.
15: comm Nag_Comm *
The NAG communication argument (see Section 3.1.1 in the Introduction to the NAG Library CL Interface).
16: fail NagError * Input/Output
The NAG error argument (see Section 7 in the Introduction to the NAG Library CL Interface).

6 Error Indicators and Warnings

NE_ALLOC_FAIL
Dynamic memory allocation failed.
See Section 3.1.2 in the Introduction to the NAG Library CL Interface for further information.
NE_BAD_PARAM
On entry, argument value had an illegal value.
NE_CONVERGENCE
The solution is not converging. See Section 9.
NE_ENUM_INT
On entry, method=Nag_Adams and iorder=2.
Constraint: if method=Nag_Adams, 3iorder6.
On entry, method=Nag_BDF and iorder=6.
Constraint: if method=Nag_BDF, 2iorder5.
NE_ENUM_INT_2
On entry, method=Nag_Adams, iorder=value and nmesh=value.
Constraint: if method=Nag_Adams, nmeshiorder-1.
On entry, method=Nag_BDF, iorder=value and nmesh=value.
Constraint: if method=Nag_BDF, nmeshiorder.
NE_INT
On entry, iorder=value.
Constraint: 2iorder6.
On entry, lwk=value.
Constraint: lwk10×nmesh+6; that is, value.
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 7.5 in the Introduction to the NAG Library CL Interface for further information.
NE_NO_LICENCE
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library CL Interface for further information.
NE_REAL
On entry, alim=value.
Constraint: alim0.0.
On entry, tol=value.
Constraint: machine precisiontol1.0.
NE_REAL_2
On entry, alim=value and tlim=value.
Constraint: tlim>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 value. If you require one more extrapolation extend the size of workspace to: lwk=value.

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 fail.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 non-stiff problems (see Section 9) you are recommended to use the Adams' method (method=Nag_Adams) of order greater than 4 (iorder>4).

8 Parallelism and Performance

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 implementation-specific information.

9 Further Comments

When solving (1), the solution of a nonlinear equation of the form
Yn-αgtn,Yn-Ψn=0, (2)
is required, where Ψn and α are constants. d05bac calls c05avc to find an interval for the zero of this equation followed by 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, d05bac will exit with fail.code= NW_OUT_OF_WORKSPACE.
In the case fail.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 gt,y in (1) with respect to its second argument is large, viz,
gt,u-gt,vLu-v. (3)
In this case, if a BDF method (method=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 (method=Nag_Adams) has been selected, you are recommended to switch to a BDF method instead.
In the case fail.code= NW_OUT_OF_WORKSPACE, then if fail.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 work[0].

10 Example

Consider the following integral equation
yt=e-t+0te-t-sys+e-ysds,  0t20 (4)
with the solution yt=lnt+e. In this example, the Adams' method of order 6 is used to solve this equation with tol=1.e−4.

10.1 Program Text

Program Text (d05bace.c)

10.2 Program Data

None.

10.3 Program Results

Program Results (d05bace.r)