NAG CL Interfaced05bec (abel1_​weak)

Settings help

CL Name Style:

1Purpose

d05bec 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.

2Specification

 #include
void  d05bec (
 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: d05bec or nag_inteq_abel1_weak.

3Description

d05bec computes the numerical solution of the weakly singular convolution Volterra–Abel integral equation of the first kind
 $f(t)+1π∫0tk(t-s) t-s g(s,y(s))ds=0, 0≤t≤T.$ (1)
Note the constant $\frac{1}{\sqrt{\pi }}$ in (1). It is assumed that the functions involved in (1) are sufficiently smooth and if
 $f(t)=tβw(t) with β>-12​ and ​w(t)​ smooth,$ (2)
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 function uses a fractional BDF linear multi-step method selected by you to generate a family of quadrature rules (see d05byc). The BDF methods available in d05bec 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 step-by-step fashion on a mesh of equispaced points. The size of the mesh is given by $T/\left(N-1\right)$, $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 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 re-evaluation of the fractional weights when d05bec is to be called several times (with the same value of $N$) within the same program with different functions.

4References

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

5Arguments

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.
userdouble *
iuserInteger *
pPointer
The type Pointer will be void *. Before calling d05bec you may allocate memory and initialize these pointers with various quantities for use by ck when called from d05bec (see Section 3.1.1 in the Introduction to the NAG Library CL Interface).
Note: ck should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d05bec. If your code inadvertently does return any NaNs or infinities, d05bec 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.
userdouble *
iuserInteger *
pPointer
The type Pointer will be void *. Before calling d05bec you may allocate memory and initialize these pointers with various quantities for use by cf when called from d05bec (see Section 3.1.1 in the Introduction to the NAG Library CL Interface).
Note: cf should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d05bec. If your code inadvertently does return any NaNs or infinities, d05bec 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.
userdouble *
iuserInteger *
pPointer
The type Pointer will be void *. Before calling d05bec you may allocate memory and initialize these pointers with various quantities for use by cg when called from d05bec (see Section 3.1.1 in the Introduction to the NAG Library CL Interface).
Note: cg should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d05bec. If your code inadvertently does return any NaNs or infinities, d05bec 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 by a previous call and are stored in rwsav.
Constraint: ${\mathbf{wtmode}}=\mathrm{Nag_InitWeights}$ or $\mathrm{Nag_ReuseWeights}$.
Note: when d05bec is re-entered with a 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: .
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: .
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×{\mathbf{iorder}}-1$, where $m\ge 1$.
9: $\mathbf{yn}\left[{\mathbf{nmesh}}\right]$double Input/Output
On entry: ${\mathbf{yn}}\left[0\right]$ must contain the value of $y\left(t\right)$ at $t=0$ (see Section 9).
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)×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 d05bec (see description of wtmode).
On exit: contains fractional weights which may be used by a subsequent call of d05bec.
11: $\mathbf{lrwsav}$Integer Input
On entry: the dimension of the array rwsav.
Constraint: ${\mathbf{lrwsav}}\ge \left(2×{\mathbf{iorder}}+6\right)×{\mathbf{nmesh}}+8×{{\mathbf{iorder}}}^{2}-16×{\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).

6Error 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.
On entry, argument $⟨\mathit{\text{value}}⟩$ 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}}=⟨\mathit{\text{value}}⟩$.
Constraint: $4\le {\mathbf{iorder}}\le 6$.
NE_INT_2
On entry, ${\mathbf{lrwsav}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{lrwsav}}\ge \left(2×{\mathbf{iorder}}+6\right)×{\mathbf{nmesh}}+8×{{\mathbf{iorder}}}^{2}-16×{\mathbf{iorder}}+1$; that is, $⟨\mathit{\text{value}}⟩$.
On entry, ${\mathbf{nmesh}}=⟨\mathit{\text{value}}⟩$ and ${\mathbf{iorder}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{nmesh}}={2}^{m}+2×{\mathbf{iorder}}-1$, for some $m$.
On entry, ${\mathbf{nmesh}}=⟨\mathit{\text{value}}⟩$ and ${\mathbf{iorder}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{nmesh}}\ge 2×{\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}}=⟨\mathit{\text{value}}⟩$.
Constraint:
On entry, ${\mathbf{tolnl}}=⟨\mathit{\text{value}}⟩$.
Constraint: .

7Accuracy

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.

8Parallelism and Performance

d05bec is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
d05bec 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 implementation-specific 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))
 $−2π k(0) g (0,y(0)) = lim t→0 f(t) t .$ (3)
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$th-order approximation can be used:
 $lim t→0 f(t) t ≃ f(hp)hp/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 $\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, d05bec computes the solution of a system of nonlinear equation for obtaining the $2p-2$ starting values. c05qdc is used for this purpose. If a failure with ${\mathbf{fail}}\mathbf{.}\mathbf{code}=$ NE_FAILED_START occurs (corresponding 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
 $Yn-αg(tn,Yn)-Ψn=0,$ (5)
is required at each step of computation, where ${\Psi }_{n}$ and $\alpha$ are constants. d05bec calls c05axc to find the root of this equation.
When a failure with ${\mathbf{fail}}\mathbf{.}\mathbf{code}=$ NE_FAILED_STEP occurs (which corresponds to an error exit from c05axc), 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{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 ill-posed or the functions in (1) are not sufficiently smooth.

10Example

We solve the following integral equations.
Example 1
The density of the probability that a Brownian motion crosses a one-sided moving boundary $a\left(t\right)$ before time $t$, satisfies the integral equation (see Hairer et al. (1988))
 $-1t exp (12-{a(t)}2/t)+∫0texp (-12{a(t)-a(s)}2/(t-s)) t-s y(s)ds=0, 0≤t≤7.$
In the case of a straight line $a\left(t\right)=1+t$, the exact solution is known to be
 $y(t)=12πt3 exp{- (1+t) 2/2t}$
Example 2
In this example we consider the equation
 $-2log(1+t+t) 1+t +∫0t y(s) t-s ds= 0, 0≤t≤ 5.$
The solution is given by $y\left(t\right)=\frac{1}{1+t}$.
In the above examples, the fourth-order BDF is used, and nmesh is set to ${2}^{6}+7$.

10.1Program Text

Program Text (d05bece.c)

None.

10.3Program Results

Program Results (d05bece.r)