nag_inteq_volterra_weights (d05bw) 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.

Syntax

[omega, lensw, sw, ifail] = d05bw(method, iorder, nomg, nwt)

[omega, lensw, sw, ifail] = nag_inteq_volterra_weights(method, iorder, nomg, nwt)

Description

nag_inteq_volterra_weights (d05bw) 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 Further Comments. For a general discussion of these methods, see Wolkenfelt (1982) for more details.

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

Parameters

Compulsory Input Parameters

1: $method$ – string (length ≥ 1)

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$ – int64int32nag_int scalar

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: $nomg$ – int64int32nag_int scalar

The number of convolution weights,

n_{w}

Constraint:

nomg \geq 1

4: $nwt$ – int64int32nag_int scalar

p

, the number of columns in the starting weights.

Constraints:

if $method ='A'$ , $nwt = iorder - 1$ ;
if $method ='B'$ , $nwt = iorder$ .

Optional Input Parameters

None.

Output Parameters

1: $omega (nomg)$ – double array: Contains the first nomg convolution weights.
2: $lensw$ – int64int32nag_int scalar: The number of rows in the weights $W_{i, j}$ .
3: $sw (ldsw, nwt)$ – double array: $ldsw = n$ .
$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 Description.
4: $ifail$ – int64int32nag_int scalar: $ifail = 0$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

Errors or warnings detected by the function:

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

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.

$ifail = - 999$: Dynamic memory allocation failed.

Accuracy

Not applicable.

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 nag_inteq_volterra_weights (d05bw). 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 nag_inteq_volterra_weights (d05bw):

(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).

Example

The following example generates the first ten convolution and thirteen starting weights generated by the fourth-order BDF method.

Open in the MATLAB editor: d05bw_example

function d05bw_example


fprintf('d05bw example results\n\n');

method = 'BDF';
iorder = int64(4);
nomg = int64(10);
nwt = int64(4);
[omega, lensw, sw, ifail] = d05bw( ...
				   method, iorder, nomg, nwt);

fprintf('\nThe convolution weights\n\n');
n = [0:double(nomg)-1]';
fprintf('%3d    %10.4f\n',[n omega]');

fprintf('\nThe weights W\n\n');
n = [1:double(lensw)]';
fprintf('%3d    %10.4f%10.4f%10.4f%10.4f\n',[n sw]');

d05bw example results


The convolution weights

  0        0.4800
  1        0.9216
  2        1.0783
  3        1.0504
  4        0.9962
  5        0.9797
  6        0.9894
  7        1.0003
  8        1.0034
  9        1.0017

The weights W

  1        0.3750    0.7917   -0.2083    0.0417
  2        0.3333    1.3333    0.3333    0.0000
  3        0.3750    1.1250    1.1250    0.3750
  4        0.4800    0.7467    1.5467    0.7467
  5        0.5499    0.5719    1.5879    0.8886
  6        0.5647    0.5829    1.5016    0.8709
  7        0.5545    0.6385    1.4514    0.8254
  8        0.5458    0.6629    1.4550    0.8098
  9        0.5449    0.6578    1.4741    0.8170
 10        0.5474    0.6471    1.4837    0.8262
 11        0.5491    0.6428    1.4831    0.8292
 12        0.5492    0.6438    1.4798    0.8279
 13        0.5488    0.6457    1.4783    0.8263

PDF version (NAG web site, 64-bit version, 64-bit version)

Chapter Contents

Chapter Introduction

NAG Toolbox