NAG Library Routine Document
D05BDF
1 Purpose
D05BDF computes the solution of a weakly singular nonlinear convolution Volterra–Abel integral equation of the second kind using a fractional Backward Differentiation Formulae (BDF) method.
2 Specification
SUBROUTINE D05BDF ( 
CK, CF, CG, INITWT, IORDER, TLIM, TOLNL, NMESH, YN, WORK, LWK, NCT, IFAIL) 
INTEGER 
IORDER, NMESH, LWK, NCT(NMESH/32+1), IFAIL 
REAL (KIND=nag_wp) 
CK, CF, CG, TLIM, TOLNL, YN(NMESH), WORK(LWK) 
CHARACTER(1) 
INITWT 
EXTERNAL 
CK, CF, CG 

3 Description
D05BDF computes the numerical solution of the weakly singular convolution Volterra–Abel integral equation of the second kind
Note the constant
$\frac{1}{\sqrt{\pi}}$ in
(1). It is assumed that the functions involved in
(1) are sufficiently smooth.
The routine uses a fractional BDF linear multistep method to generate a family of quadrature rules (see
D05BYF). The BDF methods available in D05BDF are of orders
$4$,
$5$ and
$6$ (
$\text{}=p$ say). For a description of the theoretical and practical background to these methods we refer to
Lubich (1985) and to
Baker and Derakhshan (1987) and
Hairer et al. (1988) respectively.
The algorithm is based on computing the solution
$y\left(t\right)$ in a stepbystep fashion on a mesh of equispaced points. The size of the mesh is given by
$T/\left(N1\right)$,
$N$ being the number of points at which the solution is sought. These methods require
$2p1$ (including
$y\left(0\right)$) starting values which are evaluated internally. The computation of the lag term arising from the discretization of
(1) is performed by fast Fourier transform (FFT) techniques when
$N>32+2p1$, and directly otherwise. The routine does not provide an error estimate and you are advised to check the behaviour of the solution with a different value of
$N$. An option is provided which avoids the reevaluation of the fractional weights when D05BDF is to be called several times (with the same value of
$N$) within the same program unit with different functions.
4 References
Baker C T H and Derakhshan M S (1987) FFT techniques in the numerical solution of convolution equations J. Comput. Appl. Math. 20 5–24
Hairer E, Lubich Ch and Schlichte M (1988) Fast numerical solution of weakly singular Volterra integral equations J. Comput. Appl. Math. 23 87–98
Lubich Ch (1985) Fractional linear multistep methods for Abel–Volterra integral equations of the second kind Math. Comput. 45 463–469
5 Parameters
 1: CK – REAL (KIND=nag_wp) FUNCTION, supplied by the user.External Procedure
CK must evaluate the kernel
$k\left(t\right)$ of the integral equation
(1).
The specification of
CK is:
 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 D05BDF is called. Parameters denoted as
Input must
not be changed by this procedure.
 2: CF – REAL (KIND=nag_wp) FUNCTION, supplied by the user.External Procedure
CF must evaluate the function
$f\left(t\right)$ in
(1).
The specification of
CF is:
 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 D05BDF is called. Parameters denoted as
Input must
not be changed by this procedure.
 3: CG – REAL (KIND=nag_wp) FUNCTION, supplied by the user.External Procedure
CG must evaluate the function
$g\left(s,y\left(s\right)\right)$ in
(1).
The specification of
CG is:
 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 D05BDF is called. Parameters denoted as
Input must
not be changed by this procedure.
 4: INITWT – CHARACTER(1)Input
On entry: if the fractional weights required by the method need to be calculated by the routine then set
${\mathbf{INITWT}}=\text{'I'}$ (
Initial call).
If
${\mathbf{INITWT}}=\text{'S'}$ (
Subsequent call), the routine assumes the fractional weights have been computed on a previous call and are stored in
WORK.
Constraint:
${\mathbf{INITWT}}=\text{'I'}$ or
$\text{'S'}$.
Note: when D05BDF is reentered with the value of
${\mathbf{INITWT}}=\text{'S'}$, the values of
NMESH,
IORDER and the contents of
WORK must not be changed.
 5: IORDER – INTEGERInput
On entry: $p$, the order of the BDF method to be used.
Suggested value:
${\mathbf{IORDER}}=4$.
Constraint:
$4\le {\mathbf{IORDER}}\le 6$.
 6: TLIM – REAL (KIND=nag_wp)Input
On entry: the final point of the integration interval, $T$.
Constraint:
${\mathbf{TLIM}}>10\times \mathit{machineprecision}$.
 7: TOLNL – REAL (KIND=nag_wp)Input
On entry: the accuracy required for the computation of the starting value and the solution of the nonlinear equation at each step of the computation (see
Section 8).
Suggested value:
${\mathbf{TOLNL}}=\sqrt{\epsilon}$ where $\epsilon $ is the machine precision.
Constraint:
${\mathbf{TOLNL}}>10\times \mathit{machineprecision}$.
 8: NMESH – INTEGERInput
On entry: $N$, the number of equispaced points at which the solution is sought.
Constraint:
${\mathbf{NMESH}}={2}^{m}+2\times {\mathbf{IORDER}}1$, where $m\ge 1$.
 9: YN(NMESH) – REAL (KIND=nag_wp) arrayOutput
On exit: ${\mathbf{YN}}\left(\mathit{i}\right)$ contains the approximate value of the true solution $y\left(t\right)$ at the point $t=\left(\mathit{i}1\right)\times h$, for $\mathit{i}=1,2,\dots ,{\mathbf{NMESH}}$, where $h={\mathbf{TLIM}}/\left({\mathbf{NMESH}}1\right)$.
 10: WORK(LWK) – REAL (KIND=nag_wp) arrayCommunication Array
On entry: if
${\mathbf{INITWT}}=\text{'S'}$,
WORK must contain fractional weights computed by a previous call of D05BDF (see description of
INITWT).
On exit: contains fractional weights which may be used by a subsequent call of D05BDF.
 11: LWK – INTEGERInput
On entry: the dimension of the array
WORK as declared in the (sub)program from which D05BDF is called.
Constraint:
${\mathbf{LWK}}\ge \left(2\times {\mathbf{IORDER}}+6\right)\times {\mathbf{NMESH}}+8\times {{\mathbf{IORDER}}}^{2}16\times {\mathbf{IORDER}}+1$.
 12: NCT(${\mathbf{NMESH}}/32+1$) – INTEGER arrayWorkspace
 13: IFAIL – INTEGERInput/Output

On entry:
IFAIL must be set to
$0$,
$1\text{ 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\text{ 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 $\mathbf{1}\text{ or}\mathbf{1}$ is used it is essential to test the value of IFAIL on exit.
On exit:
${\mathbf{IFAIL}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see
Section 6).
6 Error Indicators and Warnings
If on entry
${\mathbf{IFAIL}}={\mathbf{0}}$ or
${{\mathbf{1}}}$, explanatory error messages are output on the current error message unit (as defined by
X04AAF).
Errors or warnings detected by the routine:
 ${\mathbf{IFAIL}}=1$
On entry,  ${\mathbf{IORDER}}<4$ or ${\mathbf{IORDER}}>6$, 
or  ${\mathbf{TLIM}}\le 10\times \mathit{machineprecision}$, 
or  ${\mathbf{INITWT}}\ne \text{'I'}$ or $\text{'S'}$, 
or  ${\mathbf{INITWT}}=\text{'S'}$ on the first call to D05BDF, 
or  ${\mathbf{TOLNL}}\le 10\times \mathit{machineprecision}$, 
or  ${\mathbf{NMESH}}\ne {2}^{m}+2\times {\mathbf{IORDER}}1,m\ge 1$, 
or  ${\mathbf{LWK}}<\left(2\times {\mathbf{IORDER}}+6\right)\times {\mathbf{NMESH}}+8\times {{\mathbf{IORDER}}}^{2}16\times {\mathbf{IORDER}}+1$. 
 ${\mathbf{IFAIL}}=2$
The routine cannot compute the
$2p1$ starting values due to an error solving the system of nonlinear equations. Relaxing the value of
TOLNL and/or increasing the value of
NMESH may overcome this problem (see
Section 8 for further details).
 ${\mathbf{IFAIL}}=3$
The routine cannot compute the solution at a specific step due to an error in the solution of the nonlinear equation
(2). Relaxing the value of
TOLNL and/or increasing the value of
NMESH may overcome this problem (see
Section 8 for further details).
7 Accuracy
The accuracy depends on
NMESH and
TOLNL, the theoretical behaviour of the solution of the integral equation and the interval of integration. The value of
TOLNL controls the accuracy required for computing the starting values and the solution of
(2) at each step of computation. This value can affect the accuracy of the solution. However, for most problems, the value of
$\sqrt{\epsilon}$, where
$\epsilon $ is the
machine precision, should be sufficient.
In solving
(1), initially, D05BDF computes the solution of a system of nonlinear equations for obtaining the
$2p1$ starting values.
C05QDF is used for this purpose. When a failure with
${\mathbf{IFAIL}}={\mathbf{2}}$ occurs (which corresponds to an error exit from
C05QDF), you are advised to either relax the value of
TOLNL or choose a smaller step size by increasing the value of
NMESH. Once the starting values are computed successfully, the solution of a nonlinear equation of the form
is required at each step of computation, where
${\Psi}_{n}$ and
$\alpha $ are constants. D05BDF calls
C05AXF to find the root of this equation.
If a failure with
${\mathbf{IFAIL}}={\mathbf{3}}$ occurs (which corresponds to an error exit from
C05AXF), you are advised to relax the value of the
TOLNL or choose a smaller step size by increasing the value of
NMESH.
If a failure with
${\mathbf{IFAIL}}={\mathbf{2}}$ or
${\mathbf{3}}$ persists even after adjustments to
TOLNL and/or
NMESH then you should consider whether there is a more fundamental difficulty. For example, the problem is illposed or the functions in
(1) are not sufficiently smooth.
9 Example
In this example we solve the following integral equations
with the solution
$y\left(t\right)=\sqrt{t}$, and
with the solution
$y\left(t\right)=\left(1t\right)\sqrt{t}$. In the above examples, the fourthorder BDF is used, and
NMESH is set to
${2}^{6}+7$.
9.1 Program Text
Program Text (d05bdfe.f90)
9.2 Program Data
None.
9.3 Program Results
Program Results (d05bdfe.r)