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)
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
ft+1π0tkt-s t-s gs,ysds=0,  0tT. (1)
Note the constant 1π in (1). It is assumed that the functions involved in (1) are sufficiently smooth and if
ft=tβwt  with  β>-12​ and ​wt​ smooth, (2)
then the solution yt is unique and has the form yt=tβ-1/2zt, (see Lubich (1987)). It is evident from (1) that f0=0. You are required to provide the value of yt at t=0. If y0 is unknown, Section 9 gives a description of how an approximate value can be obtained.
The routine uses a fractional BDF linear multi-step 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 (=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 yt in a step-by-step fashion on a mesh of equispaced points. The size of the mesh is given by T/N-1, N being the number of points at which the solution is sought. These methods require 2p-2 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+2p-1, 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 re-evaluation 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: ck real (Kind=nag_wp) Function, supplied by the user. External Procedure
ck must evaluate the kernel kt of the integral equation (1).
The specification of ck is:
Fortran Interface
Function ck ( t)
Real (Kind=nag_wp) :: ck
Real (Kind=nag_wp), Intent (In) :: t
C Header Interface
double  ck_ (const double *t)
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 d05bef is called. Arguments denoted as Input must not be changed by this procedure.
Note: ck should not return floating-point 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: cf real (Kind=nag_wp) Function, supplied by the user. External Procedure
cf must evaluate the function ft in (1).
The specification of cf is:
Fortran Interface
Function cf ( t)
Real (Kind=nag_wp) :: cf
Real (Kind=nag_wp), Intent (In) :: t
C Header Interface
double  cf_ (const double *t)
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 d05bef is called. Arguments denoted as Input must not be changed by this procedure.
Note: cf should not return floating-point 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: cg real (Kind=nag_wp) Function, supplied by the user. External Procedure
cg must evaluate the function gs,ys in (1).
The specification of cg is:
Fortran Interface
Function cg ( s, y)
Real (Kind=nag_wp) :: cg
Real (Kind=nag_wp), Intent (In) :: s, y
C Header Interface
double  cg_ (const double *s, const double *y)
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 d05bef is called. Arguments denoted as Input must not be changed by this procedure.
Note: cg should not return floating-point 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: initwt Character(1) Input
On entry: if the fractional weights required by the method need to be calculated by the routine then set initwt='I' (Initial call).
If initwt='S' (Subsequent call), the routine assumes the fractional weights have been computed by a previous call and are stored in work.
Constraint: initwt='I' or 'S'.
Note: when d05bef is re-entered with a value of initwt='S', the values of nmesh, iorder and the contents of work must not be changed.
5: iorder Integer Input
On entry: p, the order of the BDF method to be used.
Suggested value: iorder=4.
Constraint: 4iorder6.
6: tlim Real (Kind=nag_wp) Input
On entry: the final point of the integration interval, T.
Constraint: tlim>10×machine precision.
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 9).
Suggested value: tolnl=ε where ε is the machine precision.
Constraint: tolnl>10×machine precision.
8: nmesh Integer Input
On entry: N, the number of equispaced points at which the solution is sought.
Constraint: nmesh=2m+2×iorder-1, where m1.
9: ynnmesh Real (Kind=nag_wp) array Input/Output
On entry: yn1 must contain the value of yt at t=0 (see Section 9).
On exit: yni contains the approximate value of the true solution yt at the point t=i-1×h, for i=1,2,,nmesh, where h=tlim/nmesh-1.
10: worklwk Real (Kind=nag_wp) array Communication Array
On entry: if initwt='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: lwk Integer Input
On entry: the dimension of the array work as declared in the (sub)program from which d05bef is called.
Constraint: lwk2×iorder+6×nmesh+8×iorder2-16×iorder+1.
12: nctnmesh/32+1 Integer array Workspace
13: ifail Integer Input/Output
On entry: ifail must be set to 0, -1 or 1. If you are unfamiliar with this argument you should refer to Section 4 in the Introduction to the NAG Library FL Interface 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).

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:
ifail=1
On entry, initwt=value.
Constraint: initwt='I' or 'S'.
On entry, iorder=value.
Constraint: 4iorder6.
On entry, lwk=value.
Constraint: lwk2×iorder+6×nmesh+8×iorder2-16×iorder+1; that is, value.
On entry, nmesh=value and iorder=value.
Constraint: nmesh=2m+2×iorder-1, for some m.
On entry, nmesh=value and iorder=value.
Constraint: nmesh2×iorder+1.
On entry, tlim=value.
Constraint: tlim>10×machine precision
On entry, tolnl=value.
Constraint: tolnl>10×machine precision.
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).
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).
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.
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.
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 ε, where ε 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 implementation-specific information.

9 Further Comments

Also when solving (1) the initial value y0 is required. This value may be computed from the limit relation (see Gorenflo and Pfeiffer (1991))
-2π k0 g 0,y0 = lim t0 ft t . (3)
If the value of the above limit is known then by solving the nonlinear equation (3) an approximation to y0 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 pth-order approximation can be used:
lim t0 ft t fhphp/2 . (4)
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 lim t0 ft t by analytical methods.
Also when solving (1), initially, d05bef computes the solution of a system of nonlinear equation for obtaining the 2p-2 starting values. c05qdf is used for this purpose. If a failure with ifail=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
Yn-αgtn,Yn-Ψn=0, (5)
is required at each step of computation, where Ψn and α are constants. d05bef calls c05axf to find the root of this equation.
When a failure with ifail=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 ifail=2 or 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 ill-posed 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 one-sided moving boundary at before time t, satisfies the integral equation (see Hairer et al. (1988))
-1t exp 12-at2/t+0texp -12at-as2/t-s t-s ysds=0,  0t7.  
In the case of a straight line at=1+t, the exact solution is known to be
yt=12πt3 exp- 1+t 2/2t  
Example 2
In this example we consider the equation
-2log1+t+t 1+t +0t ys t-s ds= 0,   0t 5.  
The solution is given by yt=11+t .
In the above examples, the fourth-order BDF is used, and nmesh is set to 26+7.

10.1 Program Text

Program Text (d05befe.f90)

10.2 Program Data

None.

10.3 Program Results

Program Results (d05befe.r)