NAG Library Routine Document

d05bwf (volterra_weights)


    1  Purpose
    7  Accuracy


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.


Fortran Interface
Subroutine d05bwf ( method, iorder, omega, nomg, lensw, sw, ldsw, nwt, ifail)
Integer, Intent (In):: iorder, nomg, ldsw, nwt
Integer, Intent (Inout):: ifail
Integer, Intent (Out):: lensw
Real (Kind=nag_wp), Intent (Inout):: sw(ldsw,nwt)
Real (Kind=nag_wp), Intent (Out):: omega(nomg)
Character (1), Intent (In):: method
C Header Interface
#include nagmk26.h
void  d05bwf_ (const char *method, const Integer *iorder, double omega[], const Integer *nomg, Integer *lensw, double sw[], const Integer *ldsw, const Integer *nwt, Integer *ifail, const Charlen length_method)


d05bwf computes the weights Wi,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:
0t ϕs ds h j=0 p-1 Wi,j ϕj×h + h j=p i ωi-j ϕj×h ,   0tT , (1)
with t=i×h, for i=0,1,,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 Wi,j, ωi are the starting and the convolution weights respectively. The mesh size h is determined as h=Tn, where n=nw+p-1 and nw is the chosen number of convolution weights wj, for j=1,2,,nw-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 9. For a general discussion of these methods, see Wolkenfelt (1982) for more details.


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


1:     method – Character(1)Input
On entry: the type of method to be used.
For Adams' type formulae.
For Backward Differentiation Formulae.
Constraint: method='A' or '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.
If method='A', p=iorder-1.
If method='B', p=iorder.
  • if method='A', 3iorder6;
  • if method='B', 2iorder5.
3:     omeganomg – Real (Kind=nag_wp) arrayOutput
On exit: contains the first nomg convolution weights.
4:     nomg – IntegerInput
On entry: the number of convolution weights, nw.
Constraint: nomg1.
5:     lensw – IntegerOutput
On exit: the number of rows in the weights Wi,j.
6:     swldswnwt – Real (Kind=nag_wp) arrayOutput
On exit: swij+1 contains the weights Wi,j, for i=1,2,,lensw and j=0,1,,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.
  • if method='A', ldswnomg+iorder-2;
  • if method='B', ldswnomg+iorder-1.
8:     nwt – IntegerInput
On entry: p, the number of columns in the starting weights.
  • 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 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).

Error Indicators and Warnings

If on entry ifail=0 or -1, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
On entry,method'A' or 'B'.
On entry,iorder<2 or iorder>6,
On entry,method='A' and iorder=2,
ormethod='B' and iorder=6.
On entry,method='A' and nwtiorder-1,
ormethod='B' and nwtiorder.
On entry,method='A' and ldsw<nomg+iorder-2,
ormethod='B' and ldsw<nomg+iorder-1.
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.
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.
Dynamic memory allocation failed.
See Section 3.7 in How to Use the NAG Library and its Documentation for further information.


Not applicable.

Parallelism and Performance

d05bwf is not threaded in any implementation.

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
yt=ft+0tKt,sysds,  0tT, (2)
using d05bwf. In (2), Kt,s and ft are given and the solution yt is sought on a uniform mesh of size h such that T=nh. Discretization of (2) yields
yi=fi×h+hj=0 p-1Wi,jKi,h,j,hyj+hj=piωi-jKi,h,j,hyj, (3)
where yiyi×h. We propose the following algorithm for computing yi from (3) after a call to d05bwf:
(a) Equation (3) requires starting values, yj, for j=1,2,,nwt-1, with y0=f0. These starting values can be computed by solving the linear system
yi = fi×h + h j=0 nwt-1 swij+1 K i,h,j,h yj ,   i=1,2,,nwt-1 .  
(b) Compute the inhomogeneous terms
σi = fi×h + h j= 0 nwt-1 swij+1 Ki,h,j,h yj ,   i=nwt,nwt+ 1,,n .  
(c) Start the iteration for i=nwt,nwt+1,,n to compute yi from:
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).


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

Program Text

Program Text (d05bwfe.f90)

Program Data


Program Results

Program Results (d05bwfe.r)

© The Numerical Algorithms Group Ltd, Oxford, UK. 2017