NAG CL Interface
d05bdc (abel2_weak)
1
Purpose
d05bdc 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
void 
d05bdc (
double 
(*ck)(double t,
Nag_Comm *comm),


double 
(*cf)(double t,
Nag_Comm *comm),


double 
(*cg)(double s,
double y,
Nag_Comm *comm),


Nag_WeightMode wtmode,
Integer iorder,
double tlim,
double tolnl,
Integer nmesh,
double yn[],
double rwsav[],
Integer lrwsav,
Nag_Comm *comm,
NagError *fail) 

The function may be called by the names: d05bdc or nag_inteq_abel2_weak.
3
Description
d05bdc 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 function uses a fractional BDF linear multistep method to generate a family of quadrature rules (see
d05byc). The BDF methods available in
d05bdc 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 function 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
d05bdc 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
Arguments

1:
$\mathbf{ck}$ – function, supplied by the user
External Function

ck must evaluate the kernel
$k\left(t\right)$ of the integral equation
(1).
The specification of
ck is:
double 
ck (double t,
Nag_Comm *comm)



1:
$\mathbf{t}$ – double
Input

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

2:
$\mathbf{comm}$ – Nag_Comm *
Pointer to structure of type Nag_Comm; the following members are relevant to
ck.
 user – double *
 iuser – Integer *
 p – Pointer
The type Pointer will be
void *. Before calling
d05bdc you may allocate memory and initialize these pointers with various quantities for use by
ck when called from
d05bdc (see
Section 3.1.1 in the Introduction to the NAG Library CL Interface).
Note: ck should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
d05bdc. If your code inadvertently
does return any NaNs or infinities,
d05bdc is likely to produce unexpected results.

2:
$\mathbf{cf}$ – function, supplied by the user
External Function

cf must evaluate the function
$f\left(t\right)$ in
(1).
The specification of
cf is:
double 
cf (double t,
Nag_Comm *comm)



1:
$\mathbf{t}$ – double
Input

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

2:
$\mathbf{comm}$ – Nag_Comm *
Pointer to structure of type Nag_Comm; the following members are relevant to
cf.
 user – double *
 iuser – Integer *
 p – Pointer
The type Pointer will be
void *. Before calling
d05bdc you may allocate memory and initialize these pointers with various quantities for use by
cf when called from
d05bdc (see
Section 3.1.1 in the Introduction to the NAG Library CL Interface).
Note: cf should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
d05bdc. If your code inadvertently
does return any NaNs or infinities,
d05bdc is likely to produce unexpected results.

3:
$\mathbf{cg}$ – function, supplied by the user
External Function

cg must evaluate the function
$g\left(s,y\left(s\right)\right)$ in
(1).
The specification of
cg is:
double 
cg (double s,
double y,
Nag_Comm *comm)



1:
$\mathbf{s}$ – double
Input

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

2:
$\mathbf{y}$ – double
Input

On entry: the value of the solution
$y$ at the point
s.

3:
$\mathbf{comm}$ – Nag_Comm *
Pointer to structure of type Nag_Comm; the following members are relevant to
cg.
 user – double *
 iuser – Integer *
 p – Pointer
The type Pointer will be
void *. Before calling
d05bdc you may allocate memory and initialize these pointers with various quantities for use by
cg when called from
d05bdc (see
Section 3.1.1 in the Introduction to the NAG Library CL Interface).
Note: cg should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
d05bdc. If your code inadvertently
does return any NaNs or infinities,
d05bdc is likely to produce unexpected results.

4:
$\mathbf{wtmode}$ – Nag_WeightMode
Input

On entry: if the fractional weights required by the method need to be calculated by the function then set
${\mathbf{wtmode}}=\mathrm{Nag\_InitWeights}$.
If
${\mathbf{wtmode}}=\mathrm{Nag\_ReuseWeights}$, the function assumes the fractional weights have been computed on a previous call and are stored in
rwsav.
Constraint:
${\mathbf{wtmode}}=\mathrm{Nag\_InitWeights}$ or
$\mathrm{Nag\_ReuseWeights}$.
Note: when
d05bdc is reentered with the value of
${\mathbf{wtmode}}=\mathrm{Nag\_ReuseWeights}$, the values of
nmesh,
iorder and the contents of
rwsav 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}$ – double
Input

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

7:
$\mathbf{tolnl}$ – double
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]$ – double
Output

On exit: ${\mathbf{yn}}\left[\mathit{i}1\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{rwsav}\left[{\mathbf{lrwsav}}\right]$ – double
Communication Array

On entry: if
${\mathbf{wtmode}}=\mathrm{Nag\_ReuseWeights}$,
rwsav must contain fractional weights computed by a previous call of
d05bdc (see description of
wtmode).
On exit: contains fractional weights which may be used by a subsequent call of d05bdc.

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

On entry: the dimension of the array
rwsav.
Constraint:
${\mathbf{lrwsav}}\ge \left(2\times {\mathbf{iorder}}+6\right)\times {\mathbf{nmesh}}+8\times {{\mathbf{iorder}}}^{2}16\times {\mathbf{iorder}}+1$.

12:
$\mathbf{comm}$ – Nag_Comm *

The NAG communication argument (see
Section 3.1.1 in the Introduction to the NAG Library CL Interface).

13:
$\mathbf{fail}$ – NagError *
Input/Output

The NAG error argument (see
Section 7 in the Introduction to the NAG Library CL Interface).
6
Error Indicators and Warnings
 NE_ALLOC_FAIL

Dynamic memory allocation failed.
See
Section 3.1.2 in the Introduction to the NAG Library CL Interface for further information.
 NE_BAD_PARAM

On entry, argument $\u2329\mathit{\text{value}}\u232a$ had an illegal value.
 NE_FAILED_START

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

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

On entry, ${\mathbf{iorder}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: $4\le {\mathbf{iorder}}\le 6$.
 NE_INT_2

On entry, ${\mathbf{lrwsav}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{lrwsav}}\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$.
 NE_INTERNAL_ERROR

An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact
NAG for assistance.
See
Section 7.5 in the Introduction to the NAG Library CL Interface for further information.
 NE_NO_LICENCE

Your licence key may have expired or may not have been installed correctly.
See
Section 8 in the Introduction to the NAG Library CL Interface for further information.
 NE_REAL

On entry, ${\mathbf{tlim}}=\u2329\mathit{\text{value}}\u232a$.
Constraints: ${\mathbf{tlim}}>10\times \mathit{machineprecision}$.
On entry, ${\mathbf{tolnl}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{tolnl}}>10\times \mathit{machineprecision}$.
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.
8
Parallelism and Performance
d05bdc is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
d05bdc 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 function. Please also consult the
Users' Note for your implementation for any additional implementationspecific information.
In solving
(1), initially,
d05bdc computes the solution of a system of nonlinear equations for obtaining the
$2p1$ starting values.
c05qdc is used for this purpose. When a failure with
${\mathbf{fail}}\mathbf{.}\mathbf{code}=$ NE_FAILED_START occurs (which corresponds to an error exit from
c05qdc), 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.
d05bdc calls
c05axc to find the root of this equation.
If a failure with
${\mathbf{fail}}\mathbf{.}\mathbf{code}=$ NE_FAILED_STEP occurs (which corresponds to an error exit from
c05axc), 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{fail}}\mathbf{.}\mathbf{code}=$ NE_FAILED_START or
NE_FAILED_STEP 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
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$.
10.1
Program Text
10.2
Program Data
None.
10.3
Program Results