NAG FL Interface
d05bef (abel1_weak)
1
Purpose
d05bef computes the solution of a weakly singular nonlinear convolution Volterra–Abel integral equation of the first kind using a fractional Backward Differentiation Formulae (BDF) method.
2
Specification
Fortran Interface
Subroutine d05bef ( 
ck, cf, cg, initwt, iorder, tlim, tolnl, nmesh, yn, work, lwk, nct, ifail) 
Integer, Intent (In) 
:: 
iorder, nmesh, lwk 
Integer, Intent (Inout) 
:: 
ifail 
Integer, Intent (Out) 
:: 
nct(nmesh/32+1) 
Real (Kind=nag_wp), External 
:: 
ck, cf, cg 
Real (Kind=nag_wp), Intent (In) 
:: 
tlim, tolnl 
Real (Kind=nag_wp), Intent (Inout) 
:: 
yn(nmesh), work(lwk) 
Character (1), Intent (In) 
:: 
initwt 

C Header Interface
#include <nag.h>
void 
d05bef_ ( double (NAG_CALL *ck)(const double *t), double (NAG_CALL *cf)(const double *t), double (NAG_CALL *cg)(const double *s, const double *y), const char *initwt, const Integer *iorder, const double *tlim, const double *tolnl, const Integer *nmesh, double yn[], double work[], const Integer *lwk, Integer nct[], Integer *ifail, const Charlen length_initwt) 

C++ Header Interface
#include <nag.h> extern "C" {
void 
d05bef_ ( double (NAG_CALL *ck)(const double &t), double (NAG_CALL *cf)(const double &t), double (NAG_CALL *cg)(const double &s, const double &y), const char *initwt, const Integer &iorder, const double &tlim, const double &tolnl, const Integer &nmesh, double yn[], double work[], const Integer &lwk, Integer nct[], Integer &ifail, const Charlen length_initwt) 
}

The routine may be called by the names d05bef or nagf_inteq_abel1_weak.
3
Description
d05bef computes the numerical solution of the weakly singular convolution Volterra–Abel integral equation of the first kind
Note the constant
$\frac{1}{\sqrt{\pi}}$ in
(1). It is assumed that the functions involved in
(1) are sufficiently smooth and if
then the solution
$y\left(t\right)$ is unique and has the form
$y\left(t\right)={t}^{\beta 1/2}z\left(t\right)$, (see
Lubich (1987)). It is evident from
(1) that
$f\left(0\right)=0$. You are required to provide the value of
$y\left(t\right)$ at
$t=0$. If
$y\left(0\right)$ is unknown,
Section 9 gives a description of how an approximate value can be obtained.
The routine uses a fractional BDF linear multistep method selected by you to generate a family of quadrature rules (see
d05byf). The BDF methods available in
d05bef are of orders
$4$,
$5$ and
$6$ (
$\text{}=p$ say). For a description of the theoretical and practical background related to these methods we refer to
Lubich (1987) 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
$2p2$ 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
d05bef is to be called several times (with the same value of
$N$) within the same program 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
Gorenflo R and Pfeiffer A (1991) On analysis and discretization of nonlinear Abel integral equations of first kind Acta Math. Vietnam 16 211–262
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 (1987) Fractional linear multistep methods for Abel–Volterra integral equations of the first kind IMA J. Numer. Anal 7 97–106
5
Arguments

1:
$\mathbf{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:
Fortran Interface
Real (Kind=nag_wp) 
:: 
ck 
Real (Kind=nag_wp), Intent (In) 
:: 
t 

C Header Interface
double 
ck_ (const double *t) 

C++ Header Interface
#include <nag.h> extern "C" {
double 
ck_ (const double &t) 
}


1:
$\mathbf{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
d05bef is called. Arguments denoted as
Input must
not be changed by this procedure.
Note: ck should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
d05bef. If your code inadvertently
does return any NaNs or infinities,
d05bef is likely to produce unexpected results.

2:
$\mathbf{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:
Fortran Interface
Real (Kind=nag_wp) 
:: 
cf 
Real (Kind=nag_wp), Intent (In) 
:: 
t 

C Header Interface
double 
cf_ (const double *t) 

C++ Header Interface
#include <nag.h> extern "C" {
double 
cf_ (const double &t) 
}


1:
$\mathbf{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
d05bef is called. Arguments denoted as
Input must
not be changed by this procedure.
Note: cf should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
d05bef. If your code inadvertently
does return any NaNs or infinities,
d05bef is likely to produce unexpected results.

3:
$\mathbf{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:
Fortran Interface
Real (Kind=nag_wp) 
:: 
cg 
Real (Kind=nag_wp), Intent (In) 
:: 
s, y 

C Header Interface
double 
cg_ (const double *s, const double *y) 

C++ Header Interface
#include <nag.h> extern "C" {
double 
cg_ (const double &s, const double &y) 
}


1:
$\mathbf{s}$ – Real (Kind=nag_wp)
Input

On entry: $s$, the value of the independent variable.

2:
$\mathbf{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
d05bef is called. Arguments denoted as
Input must
not be changed by this procedure.
Note: cg should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
d05bef. If your code inadvertently
does return any NaNs or infinities,
d05bef is likely to produce unexpected results.

4:
$\mathbf{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 by a previous call and are stored in
work.
Constraint:
${\mathbf{initwt}}=\text{'I'}$ or
$\text{'S'}$.
Note: when
d05bef is reentered with a value of
${\mathbf{initwt}}=\text{'S'}$, the values of
nmesh,
iorder and the contents of
work must not be changed.

5:
$\mathbf{iorder}$ – Integer
Input

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:
$\mathbf{tlim}$ – Real (Kind=nag_wp)
Input

On entry: the final point of the integration interval, $T$.
Constraint:
${\mathbf{tlim}}>10\times \mathit{machineprecision}$.

7:
$\mathbf{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 9).
Suggested value:
${\mathbf{tolnl}}=\sqrt{\epsilon}$ where $\epsilon $ is the machine precision.
Constraint:
${\mathbf{tolnl}}>10\times \mathit{machineprecision}$.

8:
$\mathbf{nmesh}$ – Integer
Input

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:
$\mathbf{yn}\left({\mathbf{nmesh}}\right)$ – Real (Kind=nag_wp) array
Input/Output

On entry:
${\mathbf{yn}}\left(1\right)$ must contain the value of
$y\left(t\right)$ at
$t=0$ (see
Section 9).
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:
$\mathbf{work}\left({\mathbf{lwk}}\right)$ – Real (Kind=nag_wp) array
Communication Array

On entry: if
${\mathbf{initwt}}=\text{'S'}$,
work must contain fractional weights computed by a previous call of
d05bef (see description of
initwt).
On exit: contains fractional weights which may be used by a subsequent call of d05bef.

11:
$\mathbf{lwk}$ – Integer
Input

On entry: the dimension of the array
work as declared in the (sub)program from which
d05bef 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:
$\mathbf{nct}\left({\mathbf{nmesh}}/32+1\right)$ – Integer array
Workspace


13:
$\mathbf{ifail}$ – Integer
Input/Output

On entry:
ifail must be set to
$0$,
$1$ or
$1$ to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of $0$ causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of $1$ means that an error message is printed while a value of $1$ means that it is not.
If halting is not appropriate, the value
$1$ or
$1$ is recommended. If message printing is undesirable, then the value
$1$ is recommended. Otherwise, the value
$0$ is recommended.
When the value $\mathbf{1}$ 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}}=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:
 ${\mathbf{ifail}}=1$

On entry, ${\mathbf{initwt}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{initwt}}=\text{'I'}$ or $\text{'S'}$.
On entry, ${\mathbf{iorder}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: $4\le {\mathbf{iorder}}\le 6$.
On entry, ${\mathbf{lwk}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{lwk}}\ge \left(2\times {\mathbf{iorder}}+6\right)\times {\mathbf{nmesh}}+8\times {{\mathbf{iorder}}}^{2}16\times {\mathbf{iorder}}+1$; that is, $\u2329\mathit{\text{value}}\u232a$.
On entry, ${\mathbf{nmesh}}=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{iorder}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{nmesh}}={2}^{m}+2\times {\mathbf{iorder}}1$, for some $m$.
On entry, ${\mathbf{nmesh}}=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{iorder}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{nmesh}}\ge 2\times {\mathbf{iorder}}+1$.
On entry, ${\mathbf{tlim}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{tlim}}>10\times \mathit{machineprecision}$
On entry, ${\mathbf{tolnl}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{tolnl}}>10\times \mathit{machineprecision}$.
 ${\mathbf{ifail}}=2$

An error occurred when trying to compute the starting values.
Relaxing the value of
tolnl and/or increasing the value of
nmesh may overcome this problem (see
Section 9 for further details).
 ${\mathbf{ifail}}=3$

An error occurred when trying to compute the solution at a specific step.
Relaxing the value of
tolnl and/or increasing the value of
nmesh may overcome this problem (see
Section 9 for further details).
 ${\mathbf{ifail}}=99$
An unexpected error has been triggered by this routine. Please
contact
NAG.
See
Section 7 in the Introduction to the NAG Library FL Interface for further information.
 ${\mathbf{ifail}}=399$
Your licence key may have expired or may not have been installed correctly.
See
Section 8 in the Introduction to the NAG Library FL Interface for further information.
 ${\mathbf{ifail}}=999$
Dynamic memory allocation failed.
See
Section 9 in the Introduction to the NAG Library FL Interface for further information.
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
(3) 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.
8
Parallelism and Performance
d05bef is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
d05bef makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the
X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the
Users' Note for your implementation for any additional implementationspecific information.
Also when solving
(1) the initial value
$y\left(0\right)$ is required. This value may be computed from the limit relation (see
Gorenflo and Pfeiffer (1991))
If the value of the above limit is known then by solving the nonlinear equation
(3) an approximation to
$y\left(0\right)$ can be computed. If the value of the above limit is not known, an approximation should be provided. Following the analysis presented in
Gorenflo and Pfeiffer (1991), the following
$p$thorder approximation can be used:
However, it must be emphasized that the approximation in
(4) may result in an amplification of the rounding errors and hence you are advised (if possible) to determine
$\underset{t\to 0}{\mathrm{lim}}}\phantom{\rule{0.25em}{0ex}}\frac{f\left(t\right)}{\sqrt{t}$
by analytical methods.
Also when solving
(1), initially,
d05bef computes the solution of a system of nonlinear equation for obtaining the
$2p2$ starting values.
c05qdf is used for this purpose. If a failure with
${\mathbf{ifail}}={\mathbf{2}}$ occurs (corresponding 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.
d05bef calls
c05axf to find the root of this equation.
When a failure with
${\mathbf{ifail}}={\mathbf{3}}$ occurs (which corresponds to an error exit from
c05axf), you are advised to either 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.
10
Example
We solve the following integral equations.
Example 1
The density of the probability that a Brownian motion crosses a onesided moving boundary
$a\left(t\right)$ before time
$t$, satisfies the integral equation (see
Hairer et al. (1988))
In the case of a straight line
$a\left(t\right)=1+t$, the exact solution is known to be
Example 2
In this example we consider the equation
The solution is given by
$y\left(t\right)=\frac{1}{1+t}$.
In the above examples, the fourthorder BDF is used, and
nmesh is set to
${2}^{6}+7$.
10.1
Program Text
10.2
Program Data
None.
10.3
Program Results