Integer, Intent (In)	::	iorder, nmesh, lwk
Integer, Intent (Inout)	::	ifail
Real (Kind=nag_wp), External	::	ck, cg, cf
Real (Kind=nag_wp), Intent (In)	::	alim, tlim, tol, thresh
Real (Kind=nag_wp), Intent (Out)	::	yn(nmesh), errest(nmesh), work(lwk)
Character (1), Intent (In)	::	method

C Header Interface

#include nagmk26.h

void

d05baf_ (
double (NAG_CALL *ck)( const double *t),
double (NAG_CALL *cg)( const double *s, const double *y),
double (NAG_CALL *cf)( const double *t),
const char *method, const Integer *iorder, const double *alim, const double *tlim, double yn[], double errest[], const Integer *nmesh, const double *tol, const double *thresh, double work[], const Integer *lwk, Integer *ifail, const Charlen length_method)

3

Description

d05baf computes the numerical solution of the nonlinear convolution Volterra integral equation of the second kind

y (t) = f (t) + \int_{a}^{t} k (t - s) g (s, y (s)) d s, a \leq t \leq T .

(1)

It is assumed that the functions involved in (1) are sufficiently smooth. The routine uses a reducible linear multi-step formula selected by you to generate a family of quadrature rules. The reducible formulae available in d05baf are the Adams–Moulton formulae of orders

3

6

, and the backward differentiation formulae (BDF) of orders

2

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,

k (t)

, 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$ – real (Kind=nag_wp) Function, supplied by the user.External Procedure

ck must evaluate the kernel

k (t)

of the integral equation (1).

The specification of ck is:

Fortran Interface

Function ck (

Real (Kind=nag_wp)	::	ck
Real (Kind=nag_wp), Intent (In)	::	t

C Header Interface

#include nagmk26.h

double

ck ( const double *t)

1: $t$ – Real (Kind=nag_wp)Input: On entry: $t$ , the value of the independent variable.

ck must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d05baf is called. Arguments denoted as Input must not be changed by this procedure.

Note: ck should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d05baf. If your code inadvertently does return any NaNs or infinities, d05baf is likely to produce unexpected results.

2: $cg$ – real (Kind=nag_wp) Function, supplied by the user.External Procedure

cg must evaluate the function

g (s, y (s))

in (1).

The specification of cg is:

Fortran Interface

Function cg (

s, y)

Real (Kind=nag_wp)	::	cg
Real (Kind=nag_wp), Intent (In)	::	s, y

C Header Interface

#include nagmk26.h

double

cg ( const double *s, const double *y)

1: $s$ – Real (Kind=nag_wp)Input: On entry: $s$ , the value of the independent variable.
2: $y$ – Real (Kind=nag_wp)Input: On entry: the value of the solution $y$ at the point s.

cg must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d05baf is called. Arguments denoted as Input must not be changed by this procedure.

Note: cg should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d05baf. If your code inadvertently does return any NaNs or infinities, d05baf is likely to produce unexpected results.

3: $cf$ – real (Kind=nag_wp) Function, supplied by the user.External Procedure

cf must evaluate the function

f (t)

in (1).

The specification of cf is:

Fortran Interface

Function cf (

Real (Kind=nag_wp)	::	cf
Real (Kind=nag_wp), Intent (In)	::	t

C Header Interface

#include nagmk26.h

double

cf ( const double *t)

1: $t$ – Real (Kind=nag_wp)Input: On entry: $t$ , the value of the independent variable.

cf must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d05baf is called. Arguments denoted as Input must not be changed by this procedure.

Note: cf should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d05baf. If your code inadvertently does return any NaNs or infinities, d05baf is likely to produce unexpected results.

4: $method$ – Character(1)Input

On entry: the type of method which you wish to employ.

$method ='A'$: For Adams' type formulae.
$method ='B'$: For backward differentiation formulae.

Constraint:

method ='A'

'B'

5: $iorder$ – IntegerInput

On entry: the order of the method to be used.

Constraints:

if $method ='A'$ , $3 \leq iorder \leq 6$ ;
if $method ='B'$ , $2 \leq iorder \leq 5$ .

6: $alim$ – Real (Kind=nag_wp)Input

On entry:

a

, the lower limit of the integration interval.

Constraint:

alim \geq 0.0

7: $tlim$ – Real (Kind=nag_wp)Input

On entry: the final point of the integration interval,

T

Constraint:

tlim > alim

8: $yn (nmesh)$ – Real (Kind=nag_wp) arrayOutput

On exit:

yn (i)

contains the most recent approximation of the true solution

y (t)

at the specified point

t = alim + i \times H

, for

i = 1, 2, \dots, nmesh

, where

H = (tlim - alim) / nmesh

9: $errest (nmesh)$ – Real (Kind=nag_wp) arrayOutput

On exit:

errest (i)

contains the most recent approximation of the relative error in the computed solution at the point

t = alim + i \times H

, for

i = 1, 2, \dots, nmesh

, where

H = (tlim - alim) / nmesh

10: $nmesh$ – IntegerInput

On entry: the number of equidistant points at which the solution is sought.

Constraints:

if $method ='A'$ , $nmesh \geq iorder - 1$ ;
if $method ='B'$ , $nmesh \geq iorder$ .

11: $tol$ – Real (Kind=nag_wp)Input

On entry: the relative accuracy required in the computed values of the solution.

Constraint:

\sqrt{ε} \leq tol \leq 1.0

, where

ε

is the machine precision.

12: $thresh$ – Real (Kind=nag_wp)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

\frac{|Y_{1} - Y_{2}|}{\max (|Y_{1}|, |Y_{2}|, |thresh|)} \leq tol,

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

thresh = 0.0

then pure relative error is measured and, if the computed solution is small and

thresh = 1.0

, absolute error is measured.

13: $work (lwk)$ – Real (Kind=nag_wp) arrayOutput

14: $lwk$ – IntegerInput

On entry: the dimension of the array work as declared in the (sub)program from which d05baf is called.

Constraint:

lwk \geq 10 \times nmesh + 6

Note: the above value of lwk is sufficient for d05baf 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

ifail = 5

6

work (1)

contains the size of lwk required for the algorithm to proceed further.

15: $ifail$ – IntegerInput/Output

On entry: ifail must be set to

0

- 1 ​ or ​ 1

. If you are unfamiliar with this argument you should refer to Section 3.4 in How to Use the NAG Library and its Documentation for details.

For environments where it might be inappropriate to halt program execution when an error is detected, the value

- 1 ​ or ​ 1

is recommended. If the output of error messages is undesirable, then the value

1

is recommended. Otherwise, if you are not familiar with this argument, the recommended value is

0

. When the value $- 1 or 1$ is used it is essential to test the value of ifail on exit.

On exit:

ifail = 0

unless the routine detects an error or a warning has been flagged (see Section 6).

6

Error Indicators and Warnings

If on entry

ifail = 0

- 1

, explanatory error messages are output on the current error message unit (as defined by x04aaf).

Errors or warnings detected by the routine:

$ifail = 1$

On entry,	$method \neq'A'$ or $'B'$ ,
or	$iorder < 2$ or $iorder > 6$ ,
or	$method ='A'$ and $iorder = 2$ ,
or	$method ='B'$ and $iorder = 6$ ,
or	$alim < 0.0$ ,
or	$tlim \leq alim$ ,
or	$tol < \sqrt{ε}$ or $tol > 1.0$ , where $ε$ is the machine precision.

$ifail = 2$

On entry,	$nmesh \leq iorder - 2$ , when $method ='A'$ ,
or	$nmesh \leq iorder - 1$ , when $method ='B'$ .

$ifail = 3$

On entry,

lwk < 10 \times nmesh + 6

$ifail = 4$: The solution of the nonlinear equation (2) (see Section 9 for further details) could not be computed by c05avf and c05azf.

$ifail = 5$: The size of the workspace lwk is too small for the required accuracy. The computation has failed in its initial phase (see Section 9 for further details).

$ifail = 6$: The size of the workspace lwk is too small for the required accuracy on the interval $[alim, tlim]$ (see Section 9 for further details).

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.
See Section 3.9 in How to Use the NAG Library and its Documentation for further information.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.
See Section 3.8 in How to Use the NAG Library and its Documentation for further information.

$ifail = - 999$: Dynamic memory allocation failed.
See Section 3.7 in How to Use the NAG Library and its Documentation for further information.

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

ifail = 5

6

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 ='A'

) of order greater than

4

(

iorder > 4

8

Parallelism and Performance

d05baf 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 routine. 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

Y_{n} - α g (t_{n}, Y_{n}) - Ψ_{n} = 0,

(2)

is required, where

Ψ_{n}

and

α

are constants. d05baf calls c05avf to find an interval for the zero of this equation followed by c05azf 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, d05baf will exit with

ifail = 5

In the case

ifail = 4

5

, you may be dealing with a ‘stiff’ equation; an equation where the Lipschitz constant

L

of the function

g (t, y)

in (1) with respect to its second argument is large, viz,

|g (t, u) - g (t, v)| \leq L |u - v| .

(3)

In this case, if a BDF method (

method ='B'

) 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 ='A'

) has been selected, you are recommended to switch to a BDF method instead.

In the case

ifail = 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 (1)

10

Example

Consider the following integral equation

y (t) = e^{- t} + \int_{0}^{t} e^{- (t - s)} [y (s) + e^{- y (s)}] d s, 0 \leq t \leq 20

(4)

with the solution

y (t) = \ln (t + 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 (d05bafe.f90)

10.2

Program Data

None.

10.3

Program Results

Program Results (d05bafe.r)

d05baf (PDF version)

D05 (inteq) Chapter Contents

D05 (inteq) Chapter Introduction

NAG Library Manual

NAG Library Routine Document

d05baf (volterra2)

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

1

Purpose

2

Specification

3

Description

4

References

5

Arguments

6

Error Indicators and Warnings

7

Accuracy

8

Parallelism and Performance

9

Further Comments

10

Example

10.1

Program Text

10.2

Program Data

10.3

Program Results