Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_inteq_volterra_weights (d05bw)

## Purpose

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 Wi,j${W}_{i,j}$ and ωi${\omega }_{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:
 t p − 1 i ∫ φ(s)ds ≃ h ∑ Wi,jφ(j × h) + h ∑ ωi − jφ(j × h),  0 ≤ t ≤ T, 0 j = 0 j = p
$∫0t ϕ(s) ds ≃h ∑ j=0 p-1 Wi,j ϕ(j×h) + h ∑ j=p i ωi-j ϕ(j×h) , 0≤t≤T ,$
(1)
with t = i × h$t=\mathit{i}×h$, for i = 0,1,,n$\mathit{i}=0,1,\dots ,\mathit{n}$, for some given constant h$h$.
In (1), h$h$ is a uniform mesh, p$\mathit{p}$ is related to the order of the method being used and Wi,j${W}_{i,j}$, ωi${\omega }_{i}$ are the starting and the convolution weights respectively. The mesh size h$h$ is determined as h = T/n$h=\frac{T}{\mathit{n}}$, where n = nw + p1$\mathit{n}={\mathit{n}}_{w}+\mathit{p}-1$ and nw${\mathit{n}}_{w}$ is the chosen number of convolution weights wj${w}_{j}$, for j = 1,2,,nw1$\mathit{j}=1,2,\dots ,{\mathit{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 [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'${\mathbf{method}}=\text{'A'}$
method = 'B'${\mathbf{method}}=\text{'B'}$
For Backward Differentiation Formulae.
Constraint: method = 'A'${\mathbf{method}}=\text{'A'}$ or 'B'$\text{'B'}$.
2:     iorder – int64int32nag_int scalar
The order of the method to be used. The number of starting weights, p$\mathit{p}$ is determined by method and iorder.
If method = 'A'${\mathbf{method}}=\text{'A'}$, p = iorder1$\mathit{p}={\mathbf{iorder}}-1$.
If method = 'B'${\mathbf{method}}=\text{'B'}$, p = iorder$\mathit{p}={\mathbf{iorder}}$.
Constraints:
• if method = 'A'${\mathbf{method}}=\text{'A'}$, 3iorder6$3\le {\mathbf{iorder}}\le 6$;
• if method = 'B'${\mathbf{method}}=\text{'B'}$, 2iorder5$2\le {\mathbf{iorder}}\le 5$.
3:     nomg – int64int32nag_int scalar
The number of convolution weights, nw${\mathit{n}}_{w}$.
Constraint: nomg1${\mathbf{nomg}}\ge 1$.
4:     nwt – int64int32nag_int scalar
p$\mathit{p}$, the number of columns in the starting weights.
Constraints:
• if method = 'A'${\mathbf{method}}=\text{'A'}$, nwt = iorder1${\mathbf{nwt}}={\mathbf{iorder}}-1$;
• if method = 'B'${\mathbf{method}}=\text{'B'}$, ${\mathbf{nwt}}={\mathbf{iorder}}$.

None.

ldsw

### 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 Wi,j${W}_{i,j}$.
3:     sw(ldsw,nwt) – double array
ldsw = n$\mathit{ldsw}=\mathit{n}$.
sw(i,j + 1)${\mathbf{sw}}\left(\mathit{i},\mathit{j}+1\right)$ contains the weights Wi,j${W}_{\mathit{i},\mathit{j}}$, for i = 1,2,,lensw$\mathit{i}=1,2,\dots ,{\mathbf{lensw}}$ and j = 0,1,,nwt1$\mathit{j}=0,1,\dots ,{\mathbf{nwt}}-1$, where n$\mathit{n}$ is as defined in Section [Description].
4:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{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${\mathbf{ifail}}=1$
 On entry, method ≠ 'A'${\mathbf{method}}\ne \text{'A'}$ or 'B'$\text{'B'}$.
ifail = 2${\mathbf{ifail}}=2$
 On entry, iorder < 2${\mathbf{iorder}}<2$ or iorder > 6${\mathbf{iorder}}>6$, or nomg < 1${\mathbf{nomg}}<1$.
ifail = 3${\mathbf{ifail}}=3$
 On entry, method = 'A'${\mathbf{method}}=\text{'A'}$ and iorder = 2${\mathbf{iorder}}=2$, or method = 'B'${\mathbf{method}}=\text{'B'}$ and iorder = 6${\mathbf{iorder}}=6$.
ifail = 4${\mathbf{ifail}}=4$
 On entry, method = 'A'${\mathbf{method}}=\text{'A'}$ and nwt ≠ iorder − 1${\mathbf{nwt}}\ne {\mathbf{iorder}}-1$, or method = 'B'${\mathbf{method}}=\text{'B'}$ and ${\mathbf{nwt}}\ne {\mathbf{iorder}}$.
ifail = 5${\mathbf{ifail}}=5$
 On entry, method = 'A'${\mathbf{method}}=\text{'A'}$ and ldsw < nomg + iorder − 2$\mathit{ldsw}<{\mathbf{nomg}}+{\mathbf{iorder}}-2$, or method = 'B'${\mathbf{method}}=\text{'B'}$ and ldsw < nomg + iorder − 1$\mathit{ldsw}<{\mathbf{nomg}}+{\mathbf{iorder}}-1$.

## Accuracy

Not applicable.

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
 t y(t) = f(t) + ∫ K(t,s)y(s)ds,  0 ≤ t ≤ T, 0
$y(t)=f(t)+∫0tK(t,s)y(s)ds, 0≤t≤T,$
(2)
using nag_inteq_volterra_weights (d05bw). In (2), K(t,s)$K\left(t,s\right)$ and f(t)$f\left(t\right)$ are given and the solution y(t)$y\left(t\right)$ is sought on a uniform mesh of size h$h$ such that T = nh$T=\mathit{n}h$. Discretization of (2) yields
 p − 1 i yi = f(i × h) + h ∑ Wi,jK(i,h,j,h)yj + h ∑ ωi − jK(i,h,j,h)yj, j = 0 j = p
$yi=f(i×h)+h∑j=0 p-1Wi,jK(i,h,j,h)yj+h∑j=piωi-jK(i,h,j,h)yj,$
(3)
where yiy(i × h)${y}_{i}\simeq y\left(i×h\right)$. We propose the following algorithm for computing yi${y}_{i}$ from (3) after a call to nag_inteq_volterra_weights (d05bw):
(a) Equation (3) requires starting values, yj${y}_{\mathit{j}}$, for j = 1,2,,nwt1$\mathit{j}=1,2,\dots ,{\mathbf{nwt}}-1$, with y0 = f(0)${y}_{0}=f\left(0\right)$. These starting values can be computed by solving the linear system
 nwt − 1 yi = f(i × h) + h ∑ sw(i,j + 1)K(i,h,j,h)yj,  i = 1,2, … ,nwt − 1. j = 0
$yi = f(i×h) + h ∑ j=0 nwt-1 swij+1 K (i,h,j,h) yj , i=1,2,…,nwt-1 .$
(b) Compute the inhomogeneous terms
 nwt − 1 σi = f(i × h) + h ∑ sw(i,j + 1)K(i,h,j,h)yj,  i = nwt,nwt + 1, … ,n. j = 0
$σi = f(i×h) + h ∑ j= 0 nwt-1 swij+1 K(i,h,j,h) yj , i=nwt,nwt+ 1,…,n .$
(c) Start the iteration for i = nwt,nwt + 1,,n$i={\mathbf{nwt}},{\mathbf{nwt}}+1,\dots ,\mathit{n}$ to compute yi${y}_{i}$ from:
 i − 1 (1 − h × omega(1)K(i,h,i,h))yi = σi + h ∑ omega( i − j + 1 )K(i,h,j,h)yj. j = nwt
$( 1 - h × omega1 K (i,h,i,h) ) y i = σ i + h ∑ 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

```function nag_inteq_volterra_weights_example
method = 'BDF';
iorder = int64(4);
nomg = int64(10);
nwt = int64(4);
[omega, lensw, sw, ifail] = nag_inteq_volterra_weights(method, iorder, nomg, nwt)
```
```

omega =

0.4800
0.9216
1.0783
1.0504
0.9962
0.9797
0.9894
1.0003
1.0034
1.0017

lensw =

13

sw =

0.3750    0.7917   -0.2083    0.0417
0.3333    1.3333    0.3333         0
0.3750    1.1250    1.1250    0.3750
0.4800    0.7467    1.5467    0.7467
0.5499    0.5719    1.5879    0.8886
0.5647    0.5829    1.5016    0.8709
0.5545    0.6385    1.4514    0.8254
0.5458    0.6629    1.4550    0.8098
0.5449    0.6578    1.4741    0.8170
0.5474    0.6471    1.4837    0.8262
0.5491    0.6428    1.4831    0.8292
0.5492    0.6438    1.4798    0.8279
0.5488    0.6457    1.4783    0.8263

ifail =

0

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

omega =

0.4800
0.9216
1.0783
1.0504
0.9962
0.9797
0.9894
1.0003
1.0034
1.0017

lensw =

13

sw =

0.3750    0.7917   -0.2083    0.0417
0.3333    1.3333    0.3333         0
0.3750    1.1250    1.1250    0.3750
0.4800    0.7467    1.5467    0.7467
0.5499    0.5719    1.5879    0.8886
0.5647    0.5829    1.5016    0.8709
0.5545    0.6385    1.4514    0.8254
0.5458    0.6629    1.4550    0.8098
0.5449    0.6578    1.4741    0.8170
0.5474    0.6471    1.4837    0.8262
0.5491    0.6428    1.4831    0.8292
0.5492    0.6438    1.4798    0.8279
0.5488    0.6457    1.4783    0.8263

ifail =

0

```