NAG Library Routine Document

D05BWF

+− Contents

1 Purpose

2 Specification

3 Description

4 References

5 Parameters

6 Error Indicators and Warnings

7 Accuracy

8 Further Comments

+− 9 Example

9.1 Program Text

9.2 Program Data

9.3 Program Results

1 Purpose

D05BWF 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.

2 Specification

SUBROUTINE D05BWF (

METHOD, IORDER, OMEGA, NOMG, LENSW, SW, LDSW, NWT, IFAIL)

INTEGER	IORDER, NOMG, LENSW, LDSW, NWT, IFAIL
REAL (KIND=nag_wp)	OMEGA(NOMG), SW(LDSW,NWT)
CHARACTER(1)	METHOD

3 Description

D05BWF computes the weights

W_{i, j}

and

ω_{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:

\int_{0}^{t} ϕ (s) d s ≃ h \sum_{j = 0}^{p - 1} W_{i, j} ϕ (j \times h) + h \sum_{j = p}^{i} ω_{i - j} ϕ (j \times h), 0 \leq t \leq T,

(1)

with

t = i \times h

, for

i = 0, 1, \dots, n

, for some given constant

h

In (1),

h

is a uniform mesh,

p

is related to the order of the method being used and

W_{i, j}

ω_{i}

are the starting and the convolution weights respectively. The mesh size

h

is determined as

h = \frac{T}{n}

, where

n = n_{w} + p - 1

and

n_{w}

is the chosen number of convolution weights

w_{j}

, for

j = 1, 2, \dots, 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 Section 8. For a general discussion of these methods, see Wolkenfelt (1982) for more details.

4 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

5 Parameters

1: METHOD – CHARACTER(1)Input

On entry: the type of method to be used.

$METHOD ='A'$: For Adams type formulae.
$METHOD ='B'$: For Backward Differentiation Formulae.

Constraint:

METHOD ='A'

'B'

2: IORDER – INTEGERInput

On entry: the order of the method to be used. The number of starting weights,

p

is determined by METHOD and IORDER.

METHOD ='A'

p = IORDER - 1

METHOD ='B'

p = IORDER

Constraints:

if $METHOD ='A'$ , $3 \leq IORDER \leq 6$ ;
if $METHOD ='B'$ , $2 \leq IORDER \leq 5$ .

3: OMEGA(NOMG) – REAL (KIND=nag_wp) arrayOutput

On exit: contains the first NOMG convolution weights.

4: NOMG – INTEGERInput

On entry: the number of convolution weights,

n_{w}

Constraint:

NOMG \geq 1

5: LENSW – INTEGEROutput

On exit: the number of rows in the weights

W_{i, j}

6: SW(LDSW,NWT) – REAL (KIND=nag_wp) arrayOutput

On exit:

SW (i, j + 1)

contains the weights

W_{i, j}

, for

i = 1, 2, \dots, LENSW

and

j = 0, 1, \dots, NWT - 1

, where

n

is as defined in Section 3.

7: LDSW – INTEGERInput

On entry: the first dimension of the array SW as declared in the (sub)program from which D05BWF is called.

Constraints:

if $METHOD ='A'$ , $LDSW \geq NOMG + IORDER - 2$ ;
if $METHOD ='B'$ , $LDSW \geq NOMG + IORDER - 1$ .

8: NWT – INTEGERInput

On entry:

p

, the number of columns in the starting weights.

Constraints:

if $METHOD ='A'$ , $NWT = IORDER - 1$ ;
if $METHOD ='B'$ , $NWT = IORDER$ .

9: IFAIL – INTEGERInput/Output

On entry: IFAIL must be set to

0

- 1 ​ or ​ 1

. If you are unfamiliar with this parameter you should refer to Section 3.3 in the Essential Introduction 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 parameter, 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'

'B'

$IFAIL = 2$

On entry,	$IORDER < 2$ or $IORDER > 6$ ,
or	$NOMG < 1$ .

$IFAIL = 3$

On entry,	$METHOD ='A'$ and $IORDER = 2$ ,
or	$METHOD ='B'$ and $IORDER = 6$ .

$IFAIL = 4$

On entry,	$METHOD ='A'$ and $NWT \neq IORDER - 1$ ,
or	$METHOD ='B'$ and $NWT \neq IORDER$ .

$IFAIL = 5$

On entry,	$METHOD ='A'$ and $LDSW < NOMG + IORDER - 2$ ,
or	$METHOD ='B'$ and $LDSW < NOMG + IORDER - 1$ .

7 Accuracy

Not applicable.

8 Further Comments

Reducible quadrature rules are most appropriate for solving Volterra integral equations (and integro-differential equations). In this section, we propose the following algorithm which you may find useful in solving a linear Volterra integral equation of the form

y (t) = f (t) + \int_{0}^{t} K (t, s) y (s) d s, 0 \leq t \leq T,

(2)

using D05BWF. In (2),

K (t, s)

and

f (t)

are given and the solution

y (t)

is sought on a uniform mesh of size

h

such that

T = n h

. Discretization of (2) yields

y_{i} = f (i \times h) + h \sum_{j = 0}^{p - 1} W_{i, j} K (i, h, j, h) y_{j} + h \sum_{j = p}^{i} ω_{i - j} K (i, h, j, h) y_{j},

(3)

where

y_{i} ≃ y (i \times h)

. We propose the following algorithm for computing

y_{i}

from (3) after a call to D05BWF:

(a)

Equation (3) requires starting values,

y_{j}

, for

j = 1, 2, \dots, NWT - 1

, with

y_{0} = f (0)

. These starting values can be computed by solving the linear system

y_{i} = f (i \times h) + h \sum_{j = 0}^{NWT - 1} SW (i, j + 1) K (i, h, j, h) y_{j}, i = 1, 2, \dots, NWT - 1 .

(b)

Compute the inhomogeneous terms

σ_{i} = f (i \times h) + h \sum_{j = 0}^{NWT - 1} SW (i, j + 1) K (i, h, j, h) y_{j}, i = NWT, NWT + 1, \dots, n .

(c)

Start the iteration for

i = NWT, NWT + 1, \dots, n

to compute

y_{i}

from:

(1 - h \times OMEGA (1) K (i, h, i, h)) y_{i} = σ_{i} + h \sum_{j = NWT}^{i - 1} OMEGA (i - j + 1) K (i, h, j, h) y_{j} .

Note that for a nonlinear integral equation, the solution of a nonlinear algebraic system is required at step (a) and a single nonlinear equation at step (c).